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Major  Department:  Mathematics 

This  dissertation  is  concerned  with  the  development  of  the  theoretical  foundation 
and  applications  of  generalized  matrix  products.  The  generalized  matrix  product,  or  p- 
product,  is  a newly  developed  matrix  operation  defined  for  heterogeneous  algebras.  This 
product  includes  the  usual  matrix  and  vector  products  of  linear  algebra  and  minimax 
algebra  as  special  cases. 

After  a brief  introduction  of  the  concept  of  a generalized  matrix  product,  several 
fundamental  properties  are  derived.  It  is  shown  that  some  properties  of  linear  trans- 
formations, such  as  the  matrix  representation  theory,  have  analogues  in  the  theory  of 
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p-products.  The  fundamental  properties  of  p-products  are  subsequently  used  to  establish 
novel  computational  methods  for  fast  transforms  used  in  digital  signal  processing. 

Specific  application  results  include  a novel  formulation  of  the  fast  Fourier  transform 
in  terms  of  the  p-product.  The  p-product  formulation  of  the  fast  Fourier  transform 
avoids  the  common  reordering  process  that  is  usually  done  either  before  or  after  a given 
computational  stage.  We  then  employ  the  p-product  as  a key  mathematical  language  in 
which  to  describe  and  analyze,  in  a unified  format,  similarities  and  differences  among 
the  Walsh,  generalized  Walsh,  and  Fourier  transforms.  Furthermore,  we  induce  a new 
algorithm  which  offers  a fast  computation  for  both  Walsh  transforms  and  does  not  require 
the  reversal  process.  Finally,  a relationship  between  the  p-product  and  the  recently 
developed  subject  wavelet  transforms  is  established.  The  p-product  formulation  provides 
for  the  first  time  a fast  version  of  the  wavelet  transform  and  also  methodologies  for 
implementation  on  parallel  computers. 
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CHAPTER  1 
INTRODUCTION 

This  dissertation  extends  the  theory  of  generalized  matrix  products  and  investigates 
the  applicability  of  generalized  matrix  operations  in  signal  processing.  The  generalized 
matrix  or  /7-product,  a newly  defined  heterogeneous  matrix  operation,  includes  the  matrix 
and  vector  products  of  linear  algebra  and  minimax  algebra  as  well  as  generalized 
convolutions.  The  dissertation  is  partitioned  as  follows:  a)  an  investigation  of  properties 
of  the  /?-product,  b)  derivation  of  a linear  transform  theory  based  on  the  /7-product, 
and  c)  derivation  of /7-product  applications  in  signal  processing.  Of  particular  interest  are 
applications  which  utilize  the  Fourier,  Walsh,  generalized  Walsh,  and  wavelet  transforms. 

The  generalized  matrix  product,  first  investigated  by  Ritter  [22],  was  defined  as  an 
operation  containing  the  matrix  product  of  linear  algebra  and  minimax  algebra  as  special 
cases.  Shortly  thereafter,  Ritter  developed  the  notion  of  the  generalized  matrix  product 
of  order  p (or  /7-product)  [21].  The  /7-product  is  a matrix  product  for  heterogeneous 
algebra  which  provides  a transformation  that  combines  the  same  or  different  values  (or 
objects)  into  values  of  a possibly  different  type  than  those  initially  used  in  the  combining 
operation.  It  represents  an  effort  to  provide  a common  mathematical  framework  in  which 
a multitude  of  ad-hoc  higher  level  computer  vision  techniques  can  be  expressed.  It  has 
been  shown  that  the  /7-product  has  certain  usual  properties  of  a matrix  product  and  can  be 
applied  to  express  various  image  processing  transforms  [24],  In  particular,  by  substituting 
specific  combining  operations  and  values  for  p,  linear  algebra  and  minimax  algebra  as  well 
as  image  algebra  operations  become  special  cases  of  this  product.  On  the  other  hand, 
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because  of  its  generality,  this  new  operation  contains  certain  unique  properties  which 
provide  some  advantages  in  parallel  processing  implementation  of  various  transforms 
[24]. 

1.1.  Background  of  a Class  of  Orthogonal  Transforms 

In  recent  years,  there  has  been  a growing  interest  in  both  the  theory  and  application 
of  orthogonal  transforms  in  the  area  of  digital  signal  processing.  This  is  primarily  due 
to  the  impact  of  high  speed  digital  computers,  rapid  advances  in  digital  technology  and 
consequent  development  of  these  advances  in  digital  signal  processors. 

One  of  the  most  fundamental  transforms  in  signal  analysis  and  processing  is  the 
Fourier  Transform  (FT).  With  rapid  advances  in  digital  computers,  the  discrete  version  of 
the  Fourier  transform  gained  in  importance.  Numerous  algorithms  have  been  developed 
for  its  efficient  computation.  Known  as  Fast  Fourier  Transforms  (FFT’s)  [27,  7,  26], 
most  of  these  algorithms  are  based  on,  or  are  variations  of,  the  Cooley-Tukey  FFT  [8], 
However,  the  foregoing  algorithms  require  two  kinds  of  operations:  namely,  the  butterfly 
and  reordering  operations.  Even  though  this  reordering  process  can  be  done  without 
arithmetic  operations,  and  involves  only  moving  a few  elements  in  storage,  it  is  still  an 
important  part  of  the  overall  cost  of  an  FFT  computation  on  most  computers.  H.  S.  Stone 
[25]  and  Cvetanovic  [9]  showed  that  these  two  kinds  of  operations  are  incompatible  with 
most  parallel  machines.  If  the  butterfly  operation  is  conflict-free,  then  the  reordering 
yields  maximum  conflict  in  the  interconnection  network,  and  vice  versa.  This  means  that 
at  least  one  of  the  two  types  of  operations  will  cause  some  problems.  The  operations 
required  for  the  reordering  can  be  up  to  O (N)  for  vectors  of  length  N.  Therefore,  the 
tedious  reordering  process  is  a major  hindrance  to  implementation  of  the  fast  Fourier 
transform  on  parallel  computers. 
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The  Walsh  Transform  (WT)  and  Generalized  Walsh  Transforms  (GWT)  are  similar 
to  the  Fourier  transform,  but  exhibit  greater  simplicity.  The  Walsh  transform  was  first 
defined  in  1923  by  J.  L.  Walsh  [28],  In  1931,  R.  E.  A.  C.  Paley  [19]  presented  a novel 
definition  of  the  Walsh  transform,  which  is  the  one  that  we  are  going  to  discuss  here. 
Paley ’s  definition  was  based  on  finite  products  of  Rademacher  functions  with  a different 
order  from  that  of  Walsh.  In  1955,  Chrestenson  extended  Paley’s  idea  and  formed 
a class  of  the  generalized  Walsh  transforms  based  on  finite  products  of  Rademacher 
functions  of  order  a,  which  takes  the  Walsh  transform  as  its  special  case  [4].  The 
Walsh  transform  can  be  computed  by  a fast  algorithm  identical  in  form  to  the  successive- 
doubling  method  given  for  the  FFT.  Due  to  the  reordering  requirement,  this  fast  method 
is  restricted  to  base-two  format  only.  It  is  probably  for  this  reason  that  fast  algorithms  for 
computing  the  generalized  Walsh  transform  have  not  been  established.  Because  of  the 
WT’s  simplicity,  efforts  to  replace  Fourier  transforms  with  Walsh  transforms  have  been 
made  in  communication  theory,  signal  processing,  image  processing,  pattern  recognition, 
etc.  [15,  2].  This  justifies  the  development  of  a general  algorithm  for  Walsh  transforms. 

It  is  well  known  that  the  Fourier  transform  decomposes  a signal  into  individual 
frequency  components  but  does  not  provide  information  as  to  when  the  frequencies 
occurred.  When  the  signal  to  be  analyzed  is  nonstationary,  a relevant  analysis  calls 
for  keeping  the  time  information  in  order  to  exhibit  its  time-varying  spectral  properties. 
The  most  straightforward  solution  involves  partitioning  of  the  signal  into  fractions  within 
which  the  stationary  assumptions  apply.  The  Gabor  Transform  (GT),  also  called  the 
Short  Time  Fourier  Transform , is  commonly  employed  in  performing  this  decomposition. 
The  GT  introduces  a time-localization  window  function , g(t  - b),  where  the  parameter  b 
translates  the  window  in  order  to  extract  local  information  from  the  Fourier  transform  of 
the  signal  over  the  entire  time  domain.  The  principal  problem  here  is  that  any  one  choice 
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of  g(t ) results  in  windows  that  are  too  wide  to  capture  all  non-stationary  behavior,  and 
too  narrow  to  capture  low-frequency  information. 

The  recently  introduced  wavelet  transform  is  an  alternative  tool  that  deals  with  non- 
stationary  signals.  Wavelet  decomposition  is  carried  out  by  means  of  a special  analysis 
function  ip,  called  the  basic  wavelet,  which  is  translated  in  time  (for  selecting  the  part  of 
the  signal  to  be  analyzed),  then  dilated  or  contracted  using  a scale  parameter  (in  order 
to  focus  on  a given  range  of  oscillations).  Wavelets  differ  from  the  Gabor  transform  in 
that  they  simultaneously  localize  a signal  and  its  Fourier  transform  with  zoom-in  and 
zoom-out  capability.  The  wavelet  transform  has  drawn  a great  deal  of  attention  from 
mathematicians  and  scientists  in  various  disciplines.  Its  numerous  applications  can  be 
found  in  [3,  6,  16,  17]. 

A well  known  fact  of  linear  algebra  is  that  linear  transforms  can  be  represented  in 
terms  of  matrix-vector  products.  However,  since  the  wavelet  transform  is  a function  of 
two  variables  (both  time  and  frequency),  it  is  expressed  with  difficulty  in  matrix  product 
form.  Even  though  Heller  et  al.  [6]  have  defined  wavelet  matrices,  they  are  used  only 
to  prove  certain  mathematical  properties.  In  order  to  be  consistent  with  other  orthogonal 
transforms  which  can  be  expressed  in  matrix  product  form  and  to  develop  more  of  its 
applications  in  digital  signal  processing,  it  is  important  to  provide  a matrix  form  of  the 
wavelet  transform. 

In  summary,  we  note  that  the  Fourier,  Walsh,  and  wavelet  transforms  exhibit  similar 
difficulties  which  can  not  be  solved  by  the  regular  matrix  product.  In  this  dissertation,  we 
show  that  the  p-product  can  be  used  to  solve  the  aforementioned  problems.  Additionally, 
we  demonstrate  the  utilizations  of  the  generalized  matrix  product  in  digital  signal 
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1.2.  Summary  of  Results 

We  outline  the  research  described  in  this  dissertation  and  the  results  that  we  have 

obtained. 

The  major  results  of  this  dissertation  are  as  follows: 

(1)  The  investigation  of  certain  matrix  product  properties  of  the  generalized  matrix 
product  and  development  of  some  of  its  unique  properties. 

(2)  Development  of  a linear  transform  theory  based  on  the  p-product. 

(3)  Derivation  and  implementation  of  a p-product  FFT  algorithm,  using  Mat  lab  3.0 
on  SUN  3 Workstations.  This  algorithm  is  designed  to  compute  the  FFT  in  a 
proper  order  without  requiring  the  reordering  process  either  after  or  before  a 
given  computational  stage. 

(4)  The  generalization  of  a new  Walsh  transform  algorithm  in  terms  of  the  p-product 
which  offers  fast  computation  for  both  Walsh  transforms  in  proper  order;  i.e., 
without  requiring  the  reverse  process. 

(5)  Evolvement  of  a class  of  p-product  matrix  transforms  over  a given  class  of  tensor 
matrix  product  transforms  [1]  based  on  derivations  of  the  p-product  Fourier  and 
Walsh  transforms. 

(6)  The  establishment  of  a relationship  between  the  p-product  and  the  wavelet 
transform  and  formulations  of  the  wavelet  transform  and  its  inverse  in  terms  of 
the  p-product.  We  provide  a simple,  fast  wavelet  transform  algorithm  that  can  be 
implemented  in  terms  of  the  p-product  on  parallel  computers. 
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In  Chapter  2,  we  define  the  generalized  matrix  product  and  its  basic  algebraic 
properties.  In  Chapter  3,  we  develop  some  linear  transformation  theories  based  on  this 
operation.  In  Chapter  4 through  Chapter  6,  we  focus  on  the  applications  of  the  p-product 
in  image  and  signal  processing.  A novel  formulation  of  the  FFT  in  terms  of  the  p- 
product  is  presented  in  Chapter  4,  which  provides  an  algorithm  that  not  only  simplifies 
the  expressions  of  the  FFT  but  also  computes  the  FFT  without  requiring  the  reordering 
process  either  after  or  before  computational  stages.  In  the  Appendix,  we  present  a 
related  computer  program.  Compared  to  the  traditional  algorithms  discussed  previously, 
only  the  butterfly  operation  is  used  in  the  p-product  FFT.  Therefore,  we  eliminate  the 
conflict  problem  on  parallel  machines,  as  well  as  the  O (N)  operations  required  for  the 
reordering  process.  Furthermore,  all  stages  of  the  p-product  FFT  are  identical,  and 
the  only  computations  required  at  each  stage  are  the  p-product  and  Hadamard  matrix 
product,  which  makes  this  algorithm  very  fast  on  either  parallel  or  sequential  machines. 
Recall  that  the  p-product  is  a key  mathematical  language  with  which  we  describe  and 
analyze,  in  a unified  format,  similarities  and  differences  among  the  Walsh,  generalized 
Walsh,  and  Fourier  transforms.  In  Chapter  5,  we  induce  a new  algorithm  which  offers  a 
fast  computation  not  only  for  both  Walsh  transforms,  but  also  eliminates  the  reordering 
process.  This  is  especially  advantageous  on  some  super  computers.  In  the  conclusion 
of  Chapter  5,  we  generate  a class  of  p-product  matrix  product  transforms  from  their 
corresponding  tensor  matrix  product  transforms.  This  demonstrates  that  any  orthogonal 
or  nonorthogonal  transforms  obtainable  from  a series  of  matrix  tensor  products  can  be 
expressed  concisely  in  terms  of  the  p-product,  and  may  be  efficiently  implemented  for 
computer  applications.  In  Chapter  6,  we  derive  a relationship  between  the  p-product 
and  the  recently  developed  wavelet  transforms.  This  provides  a simple  and  fast  wavelet 
algorithm  using  the  p-product  and  parallelism.  The  new  algorithm  decomposes  a long 


7 


summation  into  several  short  summations  and  employs  the  -product  to  carry  out  the 
computation.  When  executing  this  algorithm  on  parallel  machines,  the  computing  time 
is  g times  faster  than  the  standard  method,  where  g denotes  the  genus  of  the  wavelet 
matrix.  Conclusions  and  suggestions  for  future  research  are  given  after  Chapter  6. 


CHAPTER  2 

THE  GENERALIZED  MATRIX  PRODUCT 


In  this  chapter,  we  define  the  generalized  matrix  product  (GMP)  for  heterogeneous 
algebra  and  establish  basic  algebraic  properties  of  this  new  operation.  We  highlight  the 
close  relations  between  the  p-product  and  tensor  product,  and  demonstrate  that  the  GMP 
contains  certain  unique  properties  which  provide  some  advantages  on  parallel  computers. 
The  notation  and  properties  introduced  herein  will  be  used  throughout  this  document.  We 
define  only  those  p-product  concepts  necessary  to  describe  ideas  in  this  dissertation.  For 
a more  detailed  discussion  of  the  p-product  definition,  we  refer  the  reader  to  [21]. 


If  F is  any  set,  then  F„xm  denotes  the  set  of  all  n x m matrices  with  entries  from  F. 
Any  semigroup  (F,7)  induces,  in  a natural  way,  for  any  pair  of  positive  integers  m and 
n,  a semigroup  of  matrices  (Fnxm,  7).  The  induced  operation  on  the  matrix  semigroup 
is  defined  componentwise.  More  precisely,  if  a = (atJ)nxm  and  b = (fy)  are 
elements  of  Fnxm,  then 


The  product  37b  is  called  the  Hadamard  product. 

We  can  take  this  concept  a step  further.  If  o : E x G — ► F is  any  binary  operation, 
then  o induces  a binary  operation 


2.1.  The  Definition  of  the  Generalized  Matrix  Product 


a7b  = e = (ey) 


where 


eij  — j Vi  — 1, 2, . . . , n,j  — 1,2,...,  m. 


(1) 


(2) 
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defined  by 
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a o b = e = (eij)nxm,  where 

eij  — aij  0 bij,Vi  = 1, 2, . . . , n,j  = 1,2, ...  , m. 


(3) 


The  product  a o b is  the  Hadamard  product  for  a heterogeneous  algebra.  In  the  special 
case  where  E = G = F,  the  algebraic  structure  (F,7,o)  induces  a matrix  structure 
(F„xm,7,°)  that  behaves  very  much  like  the  structure  (F,7, o).  For  example,  if  o 
distributes  over  7 in  F,  then  o distributes  over  7 in  F„xm.  Similar  comments  can  be 
made  for  commutativity  and  associativity. 

We  will  follow  the  usual  convention  of  letting  Fm  = FiXm  and  view  Fm  as 
the  set  of  /n-dimensional  row  vectors  with  entries  from  F.  Similarly,  if  m = 1,  then 
(Fn/  = (Fix,,/  = F„xi  denotes  the  set  of  /j-dimensional  column  vectors  with  entries 


from  F. 


Let  = {1, 2, . . . , n},  p a positive  integer  dividing  both  m and  n and  m(p)  = 
y , n(p ) = j.  Define  the  following  correspondences 


cp  : Zp  x Z„+(p) 


n 


and 


by  cp(k,j ) = (k  - 1 )n(p)+j, 
where  !<_/’<  n(p),and  1 < k < p, 


rP  : Zm(p)  X Zm 


by  rp(i,k)  = ( i - l)p  + k, 
where  1 < k < p,  and  1 < i < m(p). 


(4) 


(5) 


The  Hadamard  product  defined  by  Equation  (3)  is  a componentwise  operation 
defined  only  for  matrices  of  the  same  dimension.  The  p-product  that  we  are  about 
to  discuss  is  defined  for  matrices  of  possibly  different  dimensionality  and  extends  the 
common  notions  of  vector  and  matrix  products.  In  general,  we  start  with  a heterogeneous 
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algebra  A = ({E,  F,  G},  {7,  o}),  where  o is  a binary  operation  o:ExG->F  and  (F,7) 
is  a commutative  semigroup.  For  each  quadruple  /,  m,  n,  and  q of  positive  integers  and 
each  fixed  p € Z+,  with  p dividing  both  m and  n,  we  construct  an  induced  heterogeneous 
matrix  algebra  Ap  = ({E(*m,Fjo(p)xm(p)?,  G„xm,  { ©p,7}  })  where  (Fln{p)xmfp)q,1) 
is  the  induced  matrix  semigroup  defined  earlier  and 


®p  : E/xm  x G„xg  — » ln(p)xm(p)q  (6) 

is  a binary  operation  induced  by  o and  7.  The  operation  Q)p  is  called  the  generalized 
matrix  product  of  order  p induced  by  o and  7 or  simply  the  p-product,  and  is  defined 
as  follows,  let  a = ( asj <)  e F/Xm  and  b = (&,<©  e FnXg.  Using  the  maps  rp  and  cp, 
a and  b can  be  rewritten  as 

a = {as,(i,k))lxm,  where  1 < s < /,  1 < rp(i,k)  = 3 < m,  and 

07) 

b = \b(kj)<t)nxq’  where  1 - Cp(^>i)  = *'  < n and  1 < t < q. 

The  p-product  or  generalized  matrix  product  of  a and  b is  denoted  by  a©pb,  and  is 
the  matrix 


C a©pb  GF ln^p'jxm(p)q 


(8) 


defined  by 


~ y=1  ( aa,(i,k ) 0 b{k,j),t ) 


(9) 


= K,(U)  0 o \Ptj)tt), 

where  ^ denotes  the  (s,  j)th  row  and  (i,f)th  column  entry  of  c.  Here  we  use  the 

lexicographical  order  (s,j)  < ( s',j ')  & s < s'  or  if  s = s',j  < j'.  Thus,  the  matrix  c 
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has  the  following  form: 


' ^1, 1)0.1)  • 

7i,2)0,i) 

•'  C(l.l )(!.«) 

••  7MXM) 

C(l.l)(2,l)  • 

C(1.2)(2,l) 

••  «<1.1X2, 9)  • 

••  71,2X2,?)  • 

••  71.1 )(»,«)  • 

••  7l,2)(.,0  • 

71.1) (>"(p).?) 

71.2) (m(p),?) 

c(i.n(p))(i,i)  ■ 

c(2,l)(l,l) 

••  7l.»(p)Xl.«) 

72,1)0,?) 

c(l,n(p))(2,l)  • 

C(2,l)(2,l) 

C(1,"(P))(2,?)  • 

••  72,1  )(2,?)  • 

••  7l,n<p))(i,«)  • 

• 72,1)0,0  • 

••  7i -n(p))(m(p),«) 

72,l)(m(p),?) 

C(2,n(p)Xl,l)  • 

' C(,2,n(p))(l,q) 

72,n(p))(2,l)  • 

• 72,n(p))(2,?)  • 

• 72,n(p))(i,t)  ■ 

‘ 72,"(p))(m(p),?) 

7*,i)(i,i) 

■ C(»,l)0,«) 

C(*.l)(  2,1) 

■ 7»,1)(2 ,?)  • 

• 7*,i)0,0  • 

7 *,i)(m(p),?) 

C(l.l)(l.l)  ' 

• C(,.l )(!.«) 

C('.l)(2.1) 

• C(M)(2,?)  • 

• ^UXOO  • 

7l,l)(tn(p),?) 

c(l.n(p))(l.l)  • 

• 7,.n(p))(  1.9) 

C0,n(p))(2,l)  • 

• 7,,n(p))(2>?)  • 

• 7 <,"(p))(*,‘)  • 

• 71,"(p)X”*(p),?)  ‘ 

The  entry  1)  in  the  (sj)- row  and  (i,r)-column  is  underlined  for  emphasis.  Under 
this  definition,  the  number  of  operations  required  for  computing  each  element  of  c is  p 
and  there  are  ln(p)  x m(p)q  elements  in  c.  Hence,  the  computational  complexity  of  the 
p-product  of  c = a©pb  with  a g E,xm  and  b g GnXg  is 

The  following  examples  illustrate  the  computation  of  the  GMP.  It  also  shows  that 
in  some  cases  the  p -product  can  simplify  regular  matrix  product  expressions. 


Example  2.1.1,  Let 


a = 


' a\ 

0 

02 

0 ' 

'b0' 

0 

°i 

0 

02 

and  b = 

bi 

a 3 

0 

U4 

0 

b2 

_ 0 

03 

0 

(14 

h. 

We  know  that 


xb  = 


Ol 

0 

02 

0 ‘ 

'V 

Oj&o  + 02^2 

0 

Ol 

0 

02 

v 

6t 

ai&i  + a2&3 

03 

0 

(24 

0 

A 

&2 

03  &o  + O4&2 

. 0 

03 

0 

(Z4 

63 

. 0364  -f  0463  _ 

(11) 


(12) 


However,  suppose  E = F = G = R and  7 denotes  addition  (+),  while  o denotes  multi- 
plication in  the  above  definition.  Let  a = 


«i 
a 3 


02 
a 4 


and  b be  the  same  as  above,  it 
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follows  that  / = m = 2,  n = 4,  and  g = 1.  By  taking  p=2,  one  obtains  n(p)  = 2 and 
m(p)  = 1.  According  to  the  GMP  definiton,  we  have 


a ®2^>  — 


«1 

a3 


«2 

(Z4 


®2 


' b0 ' 

aib0  -f  a2b2 

bx 

a\b\  4-  a2b$ 

b2 

a^bo  -f  a±b2 

. . 

. o.2b\  + 0463 

a x b. 


(13) 


Example  2.1.2.  Suppose  E = F = G = R,  7 denotes  addition  (+)  and  o denotes  multi- 


bi 

H 


plication.  Let  a = [01,02,03,04,05,06]  and  b = . Then  / = 1,  m = 6,  n = 2, 

and  9=1.  By  taking  p = 2,  we  obtaine  n(p)  = 1 and  m(p)  = 3.  According  to  the 
definition  of  the  GMP,  we  have 


a ©2b  — 


al  ? • a3> 


©2 


h 

b2 


(14) 


= [a\bi  + a2b2,  a2b\  -f  0462,  as&l  + aeb2]. 


Notice  that  in  the  above  examples,  n(p)  is  the  number  of  skips  on  columns  of  b 
and  m(p)  is  the  number  of  / x p blocks  in  the  matrix  a.  Therefore,  we  call  n(p)  the 
stride  of  the  GMP. 

As  mentioned  earlier,  the  /7-product  includes  the  commonly  used  matrix  and  vector 
product.  It  has  been  proved  that  those  products  can  be  obtained  by  substituting  specific 
values  for  p [21].  In  the  next  section,  we  will  discuss  some  algebraic  properties  of  the 
p-product. 

For  the  sake  of  convenience,  we  henceforth  assume  that  p is  a value  such  that  the 


various  expressions  exist. 
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2.2.  Basic  Properties  of  the  Generalized  Matrix  Product 

In  order  to  justify  the  name  the  generalized  matrix  product  and  the  notation  we  have 
used,  we  first  show  that  the  p-product  possesses  a number  of  properties  of  a product.  In 
this  section,  unless  otherwise  stated,  we  assume  (F,7,  o)  is  a ring.  We  denote  the  n x n 
identity  matrix  of  the  regular  matrix  product  by  In,  the  n x n unit  matrix  of  the  Hadamard 
matrix  product  by  ln,  and  the  general  matrix  product  generated  by  the  operations  7 and 
0 by  ©.  In  particular,  when  7 = + and  o=multiplication,  @ = x corresponds  to  the 
regular  matrix  product. 


Theorem  2.2.1.  If  a G F/Xm,  b G F nxq  and  a is  a scalar,  then 

(a)  a@p(a  o b)  = a o (a@pb), 

(b)  (a  o a)  ®pb  = a o (a©pb). 


Proof:  We  will  only  prove  (a)  since  the  proof  of  (b)  is  similar.  The  (s,  j)(i,  r)th  element 
of  a©p(aob)is 

p p 

L as,(i,k)  0 («  0 b)(kj),t  = L aSi(itk)  0 (a  o 6(Mj<) 

P 

= a 0 (Tj  0 b(k,j),t)  = the  {s,j)(i,  f)th  element  of  a 0 (a©pb). 

The  result  follows. 


(15) 


The  p-product  is  distributive  with  respect  to  7. 


Q.E.D. 


Theorem  2.2.2. 

(a)  If  a G F^xyTJ,  b G F^xm,  and  c G F uxv>  then 


(a7b)@pc  = (a©pc)7(b©pc). 


(16) 
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(b)  If  a ^ and  C6F  uxv>  then 

a®p(t>7c)  = (a®pc)7(b@pc).  (17) 

Proof:  Again  we  shall  only  prove  part  (a).  The  (s,j)(i,t) th  element  of  (a7b)  ©pc  is 

p P 

L (a76)-, (.-,*) 0 C(M.*  = F-,  (K(*.fc) 0 c(it,i),t)7(^J(,-,i)  ° <\kj),t)) 

(18) 

= the  (s,  j)(i,  t)th  element  of  (a©pc)7(b@pc). 

The  result  follows. 

Q.E.D. 

As  an  easy  corollary  we  have: 

Corollary  1.  If  a E F/Xm,  bl5  b2, . . . , b*  6 Fuxv,  and  07,  a2, . . . , a*  E F,  then 

a®P(r  («,  °b,))  =f1(a.  o(a©pb,)).  (19) 

Theorem  2.2.3.  (Associative  Law)  If  a E Fuxv,  b E Fwxl,  and  c E Fmxn>  then 

a©p(b®,c)  = (a@pb)  ©9c,  (20) 

provided  p divides  both  v and  w,  q divides  both  l and  m.  In  particular, 

ao  (b®pc)  = (ao  b)©pc,  (21) 

if  the  dimensions  of  the  matrices  are  such  that  the  various  expressions  exist. 

Proof:  Let  d = a©pb.  Then 


d G F uw(p)xv(p)h 


(22) 
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and 


P 

^(e,»/)(7,6)  =^1  (ae,(7,fc)  ° b(k,v),s)- 

Let  e = d©gc  = (a@pb)  Q)qc.  Then 


(23) 


ee  F 


uw(p)m(q)xv(p)l(p)m 


(24) 


and 


e{»,  ])(*,*)  = £j  (ds,(i,r)  ° c(r,j),t)i  (25) 

P 

where  dSj(l>)  = d(e>1|)(Ti«)  = r=i  (ae,(y,k)  0 b(k,v),s),  for  some  e,  *7,7,  6 satisfying 
•s  = (e  - l)to(p)  + rj,  and  rq(i,r ) = (7  - 1)/  + 8. 


If  f b ©?C,  then  f E F wm(q)xl(q)n  ttnd 

<7 

f(i,K)(X,p)  =I£1  (A,(A,z)  oc(x,k),^)- 

If  g = a®pf  = a®p(b®9c),  then  g € Fuu,(p)m(ff)><t)(p)/(9)fl  and 

P 

9{a,t){e,P)  =T1  {«a,(B,y)Of(yXW> 


where 

<7 

hvX),P  ~ /(t,/c)(A,/i)  — £j  bi,(X,x)  0 c(x,k), 
for  some  4,  /c,  A,  and  //  satisfying 

cp(t/,C)  = (t  - l)m(g)  + /c , and  /?  = (A  - l)n  + /i. 
In  addition,  the  (fi,f2)th  element  of  e is 

e«i<2  = e(s-l)m(?)+j,(j'-l)n+t  = e(s,j)(i,t) 

<7 

= Lj  ds,(i,r)  0 c(r,j),t 

= L CL  0 6<w) 0 


(27) 


(28) 


(29) 


(30) 


(31) 
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where 

P 

— d(e,v){l,6)  - Hi  af>(7,^)  0 hk> ?)>*’ 

with  s = (e  — l)w(p)+Ti,  and  rq(i,  r)  = (7  - l)/+<5. 
Since  rq(i,r ) = (i  — 1)<?  + r,  we  have 


8 = [*’  “ (7  “ l)/(?)  - 1]?  + r = rq(i  - (7  - 1)%),  r). 
Thus  Equation  (31)  becomes 


6<lt2  (L  fl‘.(7,*)  ° 6(*,9),(i-(7-l)/(«),r))  0 C(r,i),* 

= i£i  af>(7,fc)  ° (£j  ^(*,»?).(»-(7-l)%),r)  OC(r,j),<) 

P 

= Li  ae,(7,fc)  °/((M),i)(*‘-(7-i)/(?),0- 


(32) 


(33) 


(34) 


a = e,  C = (»7  - l)m(g)  +;,  6>  = 7,  and  /?  = [*  - (7  - l)/(?)  - l]n+*. 
Then 

p P 

=*£  <*«,(«, *)  o f(k-i)w(p)+u  = ri  ««,(»,*) 0 /(ib,c),/j 

= 9(a,C)(9,py 

However, 

9{a,C)(e,P)  — tf(a-l)w(p)m(g)+C,(0-l)/(?)ri+/? 

— 5,(f-l)w(p)m(g)+C,(7-l)/(9)n+(i-l-(7-l)/(9))n+t 


(35) 


(36) 


(37) 


^(s— l)m(g)+ji,(»— l)n+t  ~ 9titf 

Hence  = #<l(2,  and  the  result  follows. 


Q.E.D. 

If  (F,7,  o)  is  a ring  with  unity,  then  we  can  obtain  the  following  theorem. 


Theorem  2.2.4.  Assume  a e F lxm.  Then 

(a)  a ©nIn  = a = I„  ©na, 

(b)  a@p0  = 0 = 0©pa, 
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where  0 denotes  the  zero  matrix  of  appropriate  dimension. 

Proof:  Part  (b)  of  the  theorem  is  trivial.  We  only  prove  part  (a).  We  can  write  a as 
follows 

a=  [at,  a2,  am(n)]  (38) 

with  a,-  are  the  matrices  of  size  / x n.  Then 

a ® „I„  = [ai , a2, . . . , am(n)]  © „I„ 

= [at  ©In,  a2@I„,  am(n)  ®I„]  (39) 

— [ai,  a2,  am(„)]  = a. 

Similarly,  we  have 

I„©na  = a.  (40) 

Q.E.D. 

The  next  theorem  is  an  easy  corollary  of  Theorems  2.2.3  and  2.2.4. 

Corollary  1.  Let  o=multiplication,  and  7 = +.  Assume  a e F/xm,  b e FnXn.  If  b has 
an  inverse  b_1 e Fnx„  and  p=n,  then 

(a©pb)  xb_1  = a.  (41) 


Recall  from  matrix  theory  that  c = a x b is  an  upper  (or  lower)  triangular  matrix 
if  both  a and  b are  upper  (or  lower)  triangular  matrices  provided  the  expression  exists. 
Is  the  same  true  in  the  p-product?  The  following  result  provides  a necessary  condition 
for  this  property. 
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Theorem  2.2.5.  Let  a e F(xm,  b E F„xg  and  both  a and  b are  upper  triangular  matrices 
(i.e.,  a ij  = 0 for  i>j ).  If  q > n,  then  c = a©pb  is  an  upper  triangular  matrix.  Similarly 
if  both  a and  b are  lower  triangular  matrices  (i.e.,  a y = 0 for  i<j ) and  q < n,  then  c is 
a lower  triangular  matrix. 

Proof:  Assume  a G F/Xm,  b G Fnxg,  both  are  upper  triangular  matrices  and  c = a@pb. 
Then 

P 

c(sJ)(*,t)  as>(*>*)  0 b(k,j),t-  (42) 

If  both  as  (,•  >jk)  and  b(k,j),t  are  not  zeros  for  some  k,  then  C(sj)(i,t)  / 0.  It  follows  that 

s < rp(i,  k ) = (i  - l)p  + k, 

(43) 

t > cp(k,j)  = (k  - l)n(p)  + j. 

Thus 

(SJ)  = (s~  1 )n(p)+j  < ((i-l)p  + k - 1 )n(p)+j 

= (i  — l)n  + (k  — l)n(p)  + j (44) 

< (i  — l)n  + t < (i  — 1)^  + t = (i,i). 


Hence  ± 0 implies  (5,;)  < ( i,t ). 


Q.E.D. 


Let  a G 


1 0 

-1  1 


, b G [4,  2,  6,  3f,  and  suppose  7 = 4-  and  o=multiplication. 
Then  (a02b/  = [4,  2,  2,  1],  but  b'  02a'  = [4,  —2,  6,  —3].  Thus,  in  general,  the 
transpose  property  is  not  valid  for  the  p-product. 


Theorems  2.2.6-8  show  a close  relation  between  the  p-product  and  tensor  product, 
which  will  be  used  as  theoretical  building  blocks  in  later  chapters. 


Theorem  2,2.6.  Assume  (F,7,o)  is  a ring  with  unity. 

(a)  If  a e F lxm,  then 

where  k = — . 

p 

(b)  If  & (z  F mxm, 


Im  ®pa  — Ifc  ® 3, 


where  k = — . 

p 

Proof:  We  only  consider  (a),  and  (b)  can  be  proved  in  an  analogous  fashion. 
Assume  b = a®pIm,  and  c = a®  I*,  with  k = f and  Im  = (ist)mxm. 
(ti,t2)th  element  of  b is 


y 

^<!<2  = \ ) = Jlj  as,(«»  ° \x,j),tl 


for  some  fixed  i,  j,  s,  and  t.  It  follows  that 

b(s,j)(i,t)  ^ 0 = Cp(xJ)  = (x  - 1 )k  + j for  some  x <3-  x = ~~-+ i 

K 

Hence 

bhh  = j as,(«, (<-»*+!)  = a*>(«-i)p+(t-i)lfc+i  if  1 = cp(xJ),  for  s°me  x 
1 0 otherwise. 

Since  p=l  in  c = a <g>  1^, 

„ • f aa«  if  P — l, 

(“'w-r)  - - { 0 otherwise 

Let  s = a,  j — 0,  and  t = (9  — 1 )k+j  — ( i — l)m.  Then 
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(45) 


(46) 


Then  the 


(47) 


(48) 


’ (49) 


(50) 


bhh  b(s,j)(i,t)  ~ c(aJ)(9, 7)- 


(51) 
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But 


C(a./?)(0i T)  — c(a-l)k+P,(0-l)k+-f  ~ C(s-l)k+j,t+(i-l)m  ~ CU h-  (52) 


Hence  btlt2  = ctlt2.  The  result  follows. 


Q.E.D. 


Suppose  F = R,  7 = +,  and  o=multiplication,  then  corresponding  to  the 
determinant  theory  of  the  tensor  product,  we  have  the  following  result,  where  |a|  denotes 
the  determinant  of  a matrix  a. 

Theorem  2,2.7.  If  a £ Rlxm,  b £ Rnxq,  then 

(a)  a ©„b  = a x (lm(n)  ® b). 

(b)  If  in  addition,  l = m and  n — q,  then 


|a©„b|  = |a|  |br<">. 


(53) 


Proof:  We  can  write  a as 


a [ai,  a.2,  ..., 


(54) 


with  a,  are  the  matrices  of  m x n.  It  follows 


a ©nb  — [ai  x b,  a2  x b,  . . . , am(„)  x b]  = a x (diag(b,  b,  . . . , b)) 
ax  0 b). 


(55) 


For  the  tensor  product,  we  have 


|a®b|  = |a|n|b|m. 


(56) 
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Using  Equations  (55)  and  (56),  we  obtain 

|a©„b|  = |ax  (Iro(n)®b)|  = |a|  |b®Im(ri)|  = |a| |b| (57) 


Q.E.D. 


Theorem  2.2.8.  Assume  (F,  7,  o)  is  a ring  with  unity.  If  a e F/Xm  and  x is  a vector  of 
length  mn,  then  (a  ® In)x  = a©mx. 

Proof:  Since  x £ Fmnxi,  and  p = m,  the  stride  = mn/m  = n.  Thus, 

' anx\  + a\2Xn+\  + . . . + aimX(m_i)n+1 ' 

aUx2  + ai2Xn+2  + . . . + aimZ(m-l)n+2 


a\lxn  T a12x2n  T • • ■ + a\rnx  (rn—\)n+n 


= (a<g>I„)x.  (58) 


021X1  + a22X„+i  + . . . + 02m®(m-l)n+l 


_aIlxn  + o/2x2n  + . . . + a/mX(m_1)n+n 


Q.E.D. 

Additional  theorems  which  pertain  to  relationships  between  the  p-product  and 
tensor  product  will  be  discussed  in  Chapter  4.  We  next  derive  unique  properties  of 
the  generalized  matrix  product  which  have  certain  advantages  in  parallel  computing. 

Theorem  2.2,9.  //(F,7,o)  is  a ring  with  unity,  then 

(a)  Im  ©pin  =Imxn(p), 

(b)  Om  ©p0„  = Omxn(p), 


provided  p divides  both  m and  n. 
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Proof:  We  only  prove  (a).  Part  (b)  is  trivial. 

Assume  c = Im  ©pI„  = a©pb.  By  the  definition  of  the  p-product,  we  have 

P 

c(s,j)(i,t)  = £i  as,(i,k)  0 b(k,j),t- 

Thus  c(s,j)(i,t)  / 0 implies  that  there  exists  at  least  one  k,  such  that 

as,(i,k)  ± 0,  and  b(kJ)tt  ± 0. 

By  the  assumption,  we  have 

s — rp(i,  k ) = ( i - l)p  + k, 

t = cp(k,j ) = (k  - 1 )n(p)  + j. 

Recall  that 


(59) 


(60) 


(61) 


(s,j)  = {s-  1 )n{p)+j, 


(62) 


(M)  = (*  - l)n  + t. 

Substituting  Equation  (61)  into  the  above  equations,  we  obtain  (s,j)  = Hence, 


we  conclude  that 


c(s, j)(i,t)  ^ 0 =►  (s,j)  = 


(63) 


Therefore,  c = 


Q.E.D. 


Example:  Suppose  7 = +,  o=multiplication,  and  F = C.  Let 


and 


ai2 

«22 


,b  = 


bn 

b2i 


b\2 

622 


cn 

C21 


Cl  2 
C22 


(64) 


Then 

d©2c 


d - [a,  b]  = 


all  «12  b\\  b\2 

«21  a 22  621  622 


allcll  + «12C21  «11C12  + U12C22  ^llCll  + b\2C2\ 
021C11  + a22c21  a21ci2  + 0.22C22  &2lCll  + &22C21 


= [a  x c,  b x c]. 


(65) 


b\\0\2  + b\2C22 
b2lOl2  + &22C22 


This  example  leads  to  the  following  general  theorem. 


(66) 
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Theorem  2.2.10.  (Parallel  Property)  If  a € F/x„,  b <E  F/x„,  c G Fnx,, 
then 


[a,  b]  0nc  = [a  x c,  b x c]. 

In  particular,  if  ai,  a2, . . . , a*  G Fixn,  then 

[ai,a2,a3,...,afc]0nc  = [ai  x c,  a2  x c, . . . , ak  x c]. 


Proof:  The  proof  is  obvious  by  the  GMP  definition  given  in  Section  1. 


Theorem  2.2.11.  If  F = R,7  = +,  o ^multiplication,  aGF/Xfn,  and 
diagonal  matrix,  then 


a ©pd 


andi 

. . . ajpdp  aip+jdi 

• • • almdp 

a2idi 

• ••  a2pdp  a2p+idi 

• • • a2mdp 

. alldi 

. . . ajpdp  aip+idp 

• • • almdp  _ 

where 


and 


d = diag(d1,d2,...,dn). 


Proof:  The  proof  follows  directly  from  the  definition  of  the  GMP. 


Define  a function 


Ixm 


F 


Imxl 


and  F = C, 

(67) 

(68) 

Q.E.D. 
d G F nxn  O 

(69) 

(70) 

(71) 

Q.E.D. 


u : F 


(72) 


by 
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■y(a)  — [an,  ai2,  • • • ,aim,  a2i, . . . , a/m]'.  (73) 

Henceforth,  we  employ  the  notation 

co/(a)  = a(a),  and  row(a)  = (v(a))'  (74) 

to  denote  the  row  and  column  representations  of  a.  We  conclude  this  chapter  by  observing 
the  following  relationship  between  a matrix  product  and  the  /^-product. 


Theorem  2.2.12.  If  7 = +,  o=multiplication,  and  a G Ftxm,  b G Fmxg,  then 


col{  ab)  = a ©mco/(b), 

roiu(ab)  = row( a)  ®mb,  (75) 

(ab)'  = row( a)  ®TOco/(b). 


Proof:  We  only  prove  the  first  equation.  The  others  can  be  proven  in  the  same  fashion. 

Let  c = co/(b).  Then  = b ,y,  and  the  (s,j)th  element  of  a©mco/(b)  is 

given  by 


mm  m 

^ askc(k,j)  = ^ askc(k-l)q+j  = ^ ask^kj 

k= 1 k=  1 jfc=l 

= (5,  j)th  element  of  ab. 

The  result  now  follows. 


(76) 


For  the  product  of  three  matrices,  we  have 

coK a x b x c)  = (a©m(b©gco/(c))')  , 


Q.E.D. 


(77) 


qxu • 


where  c G F 


CHAPTER  3 

THE  LINEAR  ALGEBRA  INDUCED  BY  THE  GENERALIZED  MATRIX  PRODUCT 

Vectors  and  linear  transformations  play  an  important  role  in  image  processing. 
Linear  transformations  of  vector  spaces  can  be  expressed  in  terms  of  the  vector/matrix 
product.  A question  that  naturally  arises  concerns  transformations  generated  by  various 
types  of  the  generalized  matrix  product.  The  consequent  derivations  lead  to  a definition  of 
p-product  transformations.  Under  such  notation,  we  establish  the  matrix  representation  of 
the  p-product  transformation,  which  includes  the  representation  of  a linear  transformation 
as  its  special  case.  In  the  final  section,  we  develop  a new  vector  space  called  the  p-product 
vector  space. 

Since  vector  spaces  and  linear  transformations  mainly  deal  with  additions  and  scalar 
multiplications,  we  denote  7 by  + and  o by  multiplication  throughout  the  remainder  of 
this  chapter. 

3.1.  p-product  Transformations 

In  this  section,  after  briefly  introducing  the  concepts  of  vector  spaces,  we  define 
the  notion  of  p-product  transformation. 

Definition  3.1.1.  A vector  space  V over  the  field  F,  denoted  by  V(F),  is  an  additive 
abelian  group  V together  with  an  operation  called  scalar  multiplication  of  each  element 
of  V by  each  element  of  F on  the  left,  such  that  Va,  /?  e F and  v,  w e V,  the  following 
five  conditions  are  satisfied: 

a)  a • v e V 
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b)  a • (/?  • v)  = (a  • /?)  • v 

c)  (a  + /?)  • v = (a  • v)  + (/3  • v) 

d)  a • (v  + w)  = (a  • v)  + (a  • w) 

e)  1 • v = v. 

The  elements  of  V are  called  vectors  and  the  elements  of  F are  called  scalars. 

If  the  field  F is  clear  from  the  context  of  discussion,  then  it  is  customary  to  use 
the  symbol  V instead  of  V(F). 

It  is  obvious  that  a set  (Fmxn,  +,  •)  of  all  m x n matrices  with  entries  from  F is 
a vector  space. 

In  the  study  of  vector  spaces,  the  most  important  types  of  mappings  are  linear 
transformations. 

Definition  3.1.2,  A linear  transformation  or  linear  operation  of  a vector  space  V(F) 
into  a vector  space  W(F)  is  a function  L : V(F)  ->  W(F)  which  satisfies 

L(a  • u + 0 • v)  = a • L(u)  + f3  • L(v),  (78) 

for  all  u,v  € V(F)  and  for  all  scalars  a,  ft  6 F. 

If  a is  any  m x n matrix  with  entries  from  F,  we  define  a linear  operator 
La  : F"  ->  Fm  by 

La(x)  = a x x,  (79) 

for  each  x £ F".  The  operator  La  is  a linear  transformation  since 
La(a  - x + /3-y)  = ax(a-x  + /?-y) 

= a x (a  ■ x)  + a x (0  ■ y)  = a • (a  x x)  + 0 ■ (a  x y)  (80) 

= a • La(x)  + 0 ■ La(y). 

We  call  La  a left-multiplication  transformation. 

Similar  to  this,  we  give  the  following  definition  for  a left  p-product  transformation. 
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Definition  3.1,3.  Let  a E F/Xm.  A left  p-product  transformation  of  the  vector  space 
(F»xg,  +,  ■)  into  a vector  space  (Fuxt),  +,  •)  is  a function  ga  : Fnxg  -♦  Fuxt;  defined  by 

ga(b)  = a©pb,  (81) 

for  all  b E F nxq,  where  u = l x n(p)  and  v = m{p ) x q. 

Note  that  the  left  p-product  transformation  is  a linear  transformation.  In  fact,  it 
follows  from  Corollary  1 of  Theorem  2.2.2  that  for  all  a,/3  e F, 

ga(ob  + /?c)  = a ©p(ab  + /3c)  = a ©p(ab)  + a ®p(/?c) 

(82) 

= a(a  ©pb)  + /?(a  ©pc)  = oga(b)  + ^9ga(c). 

Hence,  all  of  the  important  properties  of  linear  transformations  are  valid  for  the  left 
p-product  transformations. 

Corresponding  to  left  p-product  transformations,  a definition  of  a right  p-product 
transformation  is  given  below. 

Definition  3.1.4.  Let  a eF;xm.  A right  p-product  transformation  of  a vector  space 
F„xg  into  a vector  space  Fux„  is  a function  fa  : Fnxg  -+  Fuxv  defined  by 

fa(b)  = b ©pa,  (83) 

for  all  b E Fnxg,  where  u = / x n(p),  and  v = m(p)  x q. 

It  is  clear  that  right  p-product  transformations  have  properties  which  are  similar  to 
those  of  left  p-product  transformations. 
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3.2.  The  Matrix  Representation  of  p-product  Transformations 

We  have  given  the  definition  of  the  p-product  transformation.  We  shall  now  embark 
upon  one  of  the  most  useful  approaches  to  the  analysis  of  a p-product  transformation  on  a 
finite-dimensional  vector  space,  the  matrix  representation  of  a p-product.  In  fact,  similar 
to  linear  transformations,  we  develop  a one-to-one  correspondence  between  matrices 
and  p-product  transformations  that  will  allow  us  to  utilize  properties  of  one  to  study 
properties  of  the  other. 

In  the  following  discussion,  assume  F is  a ring  with  unity  and  Fnx?  and  Fuxu 
are  the  vector  spaces  over  F.  In  the  last  section  we  have  shown  that  for  each 
a6Fixra,  there  exists  a left  p-product  transformation  ga  from  Fnx,  into  Fux„,  with 
« = l x n(p)  and  v = m(p)  x q,  such  that  ga(b)  = a©pb  for  every  b £ FnXq.  It  is 
easy  to  verify  that  ga  is  a homomorphism. 


Theorem  3.2,1,  Let  u,  v,  l,  m,  n,  q,  and  p be  positive  integers  with  u = / x 
n(p)  and  v = m(p)  x q.  The  set  Homp(FnX9,  Fuxt))  = {ga  : a £ F/xm}  of  homomor- 
phisms  from  Fnx?  to  Fuxv  is  an  abelian  group  under  addition. 

Proof.  Assume  a,  b,  c £ F/Xm  and  let  ga,  g^,  gc  £ Homp(F nxq,  F uxv).  Then  ga  + g^  £ 
Homp(FnX(?,Fuxw)  since 

(ga  + gb)(x)  = ga(x)  + gb(x)  = a ©px  + b ©px 


= (a  + b)  ©px  = ga+b(x), 

for  every  x £ Fnxg.  Thus  ga  + gb  = ga+b-  But  ga+b  £ Homp(F„X9, 
F/Xr„  is  an  abelian  group  under  addition.  Clearly, 


(84) 

Fux„)  since 


ga  + (gb  + gc)  = (ga  + gb)  + gc, 


(85) 


since 
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ga  + (gb  + gc)  - ga+(b-fc)  = g(a+b)+c  = (ga  + gb)  + gc- 


Similarly, 


ga  "1  gb  gb  ga? 


(86) 


(87) 


since 


ga+b  — gb+a- 


(88) 


Finally,  it  is  easy  to  verify  that  the  identity  element  is  0 =g0  and  the  inverse  of  ga  is 
g_a.  Hence  Homp(Fnxg,  Fux„)  is  an  abelian  group  under  addition. 


Q.E.D. 


Note  that  for  Va,  /?  G F andga,gb  G Homp(Fnxg,  Fox„),  we  have 
oga  G Homp(F nx?,Fux„), a(/?ga)  = (a(3) ga, 


and  a(ga  + gb)  = aga  + agb. 

Hence,  it  is  not  difficult  to  prove  that  Homp(F„X9,  Fux„)  is  a vector  space. 

The  next  theorem  provides  for  an  important  relationship  between  Homp(F„x?, 
and  the  set  of  matrices  F/Xm. 


(89) 


F uxt>) 


Theorem  3.2.2.  Let  F be  a ring  with  identity.  Let  F/xm,FnX9,  and  Fux„  be  vector 
spaces  over  F with  u = l x n(p)  and  v = m(p)  x q for  some  fixed  integer  p dividing 
both  m and  n.  Then  there  exists  a vector  space  isomorphism 


i()  . Homp(F nxq,  Fux„ 


x m • 


(90) 
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Proof:  Let  ip  : F/Xm  — ► Homp(Fnxg,  Fux„)  be  defined  by 


V’(a)  = ga. 


(91) 


The  function  ip  is  well-defined  since 


ip(a)  ± V’(b)  =>  ga  ± gb  =>  a ©px  ± b ©px,  for  some  x G F„xg, 

=>  (a  - b)  ©px  ± 0 =4*  a / b. 

Also,  ip  preserves  the  operation  of  addition  since 


(92) 


ip(a  + b)  = ga+5  = ga  + g5  = ip(a)  + ip( b). 


(93) 


It  is  easy  to  verify  that  ip  is  onto  by  using  the  definition  of  Homp(F„xg,  Fux„). 

We  show  that  ip  is  one-to-one.  If  ip  (a)  = 0,  then  ga  = 0.  Thus  ga(b)  = 0 for  every 
b G Fnxg;  i.e.  a ©pb  = 0 for  every  b.  Assume  a^O;  i.e.  a ^ 0 for  some  i,j.  Let 
c = a ©pb.  Then  by  the  definition  of  the  ^-product,  we  have 


If  &s,(i,k)  7^  0 f°r  some  fixed  s,  i and  k,  then  we  can  choose  a matrix  b,  such  that 
b(jtj),t  7^  0 for  fixed  j,  t.  Obviously  b G Fnxg.  Hence  ^ ^ 0,  which  is  a 

contradiction.  Thus  a=0.  Therefore,  ip  is  an  isomorphism. 

Q.E.D. 

If  V is  an  n dimensional  vector  space  over  F and  W is  an  m dimensional  vector 
space  over  F,  then  we  have  the  following  corollary. 

Corollary  L Suppose  pm  = In  and  for  a G F/Xp,  let  ga  : V — > W,  be  defined  by 


p 


(94) 


jt=i 


ga(v)  = a©pV. 


(95) 
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//■Homp(V,W)  = {ga  : a £ F/Xp,andp  = ^},  then  the  function 

V>  : {FZxp,+}  - {Homp(V , W),  +},  (96) 

defined  by 

V>(a)  = ga,  (97) 

is  an  isomorphism. 

We  call  F/Xp  the  representation  of  Homp(V,  W),  where  / = Notice  that  dif- 
ferent p give  different  representations  for  Homp(V,  W).  For  example,  if  p-n,  then  l=m. 
Thus  Hom„(V,  W)  ss  Fmxn  under  the  normal  matrix  product,  which  is  the  representation 
of  linear  transformations  from  V to  W.  If  p= 1,  then  Homi(V,  W)  « F/xl  = F-Xl  under 

n 

the  tensor  product.  The  question  that  has  been  raised  here  is  that  since  in  the  theory  of 
linear  transformations,  each  homomorphism  can  be  represented  by  different  matrices 
relative  to  different  bases,  does  the  same  property  hold  for  the  theory  of  p-product 
transformations?  We  plan  to  determine  to  answer  this  question  in  a future  research. 

We  conclude  this  section  with  the  following  composition  theorem. 

Theorem  3.2,3.  Let  U be  a q dimensional  vector  space  over  F with  q < n,  and  let  V 
and  W be  the  same  as  above.  If  g £ Homp(V,  W)  corresponds  to  the  l x p matrix  a 
with  p = and  f £ Homp(U,  V)  corresponds  to  the  s x p matrix  b with  s = then 
g o f £ Homp(U,  W)  corresponds  to  the  y x ^ matrix  a ©pb. 

Proof:  Let  u £ U.  Then 

gof(u)  = g(f(u))  = g(b©pu)  = a©p(b©pu)  = (a ©pb)  ©pu.  (98) 


Q.E.D. 


3.3.  Generalized  Vector  Space  Products 


32 


Recall  that  the  tensor  product  V®W  of  two  vector  spaces  V and  W of 
dimensions  n and  m,  respectively,  is  a vector  space  with  dimension  mn.  This  raises 
a question  as  to  the  possibility  of  a p-product  of  V and  W,  which  is  discussed  in  this 
section. 


Example  1:  Let  V = R2  with  elementary  basis  {ei,e2},  W = R4  with  elementary  basis 
{ei,  e2,  e3,  64},  and  p = 2 with 


©2  ; v x W'  ->  V', 


(99) 


defined  by 


(x,y)  = x ®2y. 


(100) 


Then  we  have 


and 


[1  0]  ©2 


[1  0]  ©2 


[0  1]  ©2 


T 

'o' 

0 

0 

= 

r 

0 

5 [1  0]  ©2 

1 

0 

0_ 

_0_ 

'o' 

'o' 

0 

'o' 

o 

1 

— 

0 

5 [1  0]  ©2 

0 

0 

1_ 

V 

'o' 

0 

'o' 

1 

0 

0 

> [0  1]  ©2 

0 

0_ 

0 

'o' 

'o' 

0 

r 

o 

1 

= 

0 

> [0  1]  ©2 

0 

0_ 

_1 

0 

0 


0 

1 


(101) 


(102) 


[0  1]  ©2 
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Note  that  if  we  denote  V ©2W  = V',  then 

M 


(103) 


0 


and 


0 


(104) 


0 


form  a basis  of  V @2W.  Thus  V ©2W  is  a vector  space  of  dimensions  2 with  a basis 
generated  by  some  elements  of  the  bases  V and  W.  The  next  theorem  generalizes  this 
concept. 

Let  V be  an  ^-dimensional  vector  space  over  F with  basis  {vi,  v2, . . . , v„}  and 
W be  an  m-dimensional  vector  space  over  F with  basis  {wi,w2, . . . wm}. 

Theorem  3.3.1.  If  p divides  both  m and  n,  then  V ©pW  is  a vector  space  of  dimension 
m(p)  x n(p)  with  basis  given  by 


v(n(p)-l)*p+l  ©pWi,  v(n(p)_1)tp+1©pw2)  ...  v(n(p)_1)+p+1  ©pwm(p)}. 

Proof:  Assume  that  for  1 < i < n(p),  1 < j < m{p),  there  exist  a, j 6 F such  that 
atJ  f 0 for  some  i and  j,  and 

an  (vi  ©pWl)  + a2i  (vi  ©pw2)  -(-...  + am(p)1  (vi  ©pwm(p))  + 
ai2(vp+i  ©pWl)  + . . . -I-  am(p)2(vp+i  ©pwm(p))  + . . .+ 

aln(p)  (v(n(p)  — l)*p+l  ©pwl)  T am(p)n(p)  (v(n(p)— l)*p+l  ©pwm(p)) 

= (anvi  + ai2vp+i  -I-  . . . + ai„(p)V(n(p)_1)+p+1)  ©pwi  + 


{vi  ©pwi, 

vP+l  ©pwl 


Vl  ©pW2, 
VP+1  ©pw2 


Vl  ©pWm(p), 

Vp+1 


(105) 


(a21V!  + a22vp+1  + . . • + a2n(p)V(„(p)_i)+p+i)  ©pw2  + . . . + 

(am(p)lvl  + . . . + am(p)n(p)v(n(p)  — i)*p-j_i)  ©pwm(p)  = 0- 


(106) 
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Write  w,  as  a matrix  (w‘t)pxm(p)  for  each  i.  That  is,  if  w,  = (c1?  c2, . . . , cm)',  then 


«t) 


pxm(p) 


Cl  C2 

Cm(p)+ 1 cm(p)+ 2 


cm(p) 

c2m(jp) 


(107) 


Let 


- c(p-l)m(p)+l  ‘ • • • cpm(p ) J 

Since  { wi , w2 , . . . , wm  } is  linearly  independent,  so  is  j (w*t) , (w^) , . . 

anvi  + ai2vp+1  . . . + ai„(p)V(„(p)_1)*p+1 

021V1  + a22Vp+i  + . . . + a2n(p)v(n(p)-l)*p+l 

G = 

.am(p)  1V1  + • • • + am(p)n(p)v(n(p)— l)*p+l 

Then  G is  an  m(p)  x n matrix.  We  rewrite  G as 


,(w”W)}. 


(108) 


G = 


511 

512 

521 

522 

- 5m(p)  1 

5m(p)2 

9ln(p) 

92n(p) 


(109) 


• • 9m(p)n(p)  - 

where  each  gij  is  a matrix  of  size  1 x p.  Using  Theorem  2.2.10,  we  may  write 
Equation  (106)  as 


[911,912,  ■ ■ • ,5ln(p)]  ©p(w^)  + 

[<721,  <722,  • • • , <72n(p)]  ©p(wst)  + •••  + 
[9m(p)l,  9m(p)2,  • • ■ , 9m(p)n(p)\  ©p  ^ = 

[511  (wjt)  512  (wjt)  . . . 5l«(p)  (wat)]  + • • • + 

Mp)i(W™(P))  ••■Mp)n(p)(W™(P))j  = 
[911  (wi)  + . . . + gm{p)1  (w-(p))  + 

"h  • • • 5ln(p)(wjt)  + • • ■ + 9m(p)n(p)  st  )] 
= 0. 
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This  implies  that  gT^  = 0 for  all  r and  k.  Hence  = 0 for  all  i and  j since 
{vi,  V2, . . . , v„}  is  linearly  independent,  which  is  a contradiction. 

Q.E.D. 

Note  that  if  we  take  p = 1,  then  V ©pW  is  a vector  space  of  dimensions  mn, 
with  a basis 

{Vi  ® W!,Vi  0 w2,...,vi  (8)  wm> 

V2  0W1,V2®W2 V2®W  m, 

(111) 


v„  ® W!,V„  ® W2,...V„  0 wm}, 

which  is  exactly  the  tensor  product  space.  Therefore,  by  substituting  specific  values  of 
p,  we  can  obtain  different  p-product  spaces  which  contain  the  tensor  product  space  as 
a special  case. 

For  any  element  x E V and  y e W,  we  have 


n m 

X = ^2  a,v,-,  and  y = bjWj,  (112) 

t=l  j= l 

with  a,-,  bj  E F.  Then 

n m 

x a*v«  ®p 

i=l  j=l 

n m 

= a'b3  (v*  ©Pwi)  (1 13) 

i=l  j= 1 

= SC«j(V<©pWj)- 


In  order  to  be  consistent  with  the  tensor  product  notation,  we  write  each  element 
z E V ©pW  as 

k 

a»(x,  ©py,),  with  a,-  E F,  x,  E V,  and  y,  E W. 

i=l 


(114) 
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Definition  3.3.2.  An  element  z e V ©p  W is  said  to  be  of  rank  k if  z = J2  (x,  ® y,-) , 

»'=i 

where  {xi, x2, . . . , x*}  and  {yi, y2, . . . , y^}  are  linearly  independent  and  k denotes  the 
minimum  number  of  terms  of  the  form  x ©py  in  which  z can  be  expressed. 

The  next  theorem  provides  another  view  of  V ©pW. 


Theorem  3.3.3.  If  Fm(p)xrj(p)  denotes  the  vector  space  ofm(p)  x n(p)  matrices  over  F, 
then  V ®pW  is  isomorphic  to  Fm(p)xn(p). 

Proof:  Define  a map  <p  by 


E a 


^(Vi©?Wj)  = 


if  i = 1, 
Ej(i-)-i)  if  i = p + l, 


,EMp)  if  i = {n(p) -l)p  + l, 

where  EJt  is  the  matrix  with  1 in  the  (j,i)  position  and  0 elsewhere, 
the  basis  of  V @pW  one  to  one  and  onto  the  basis  of  Fra(p)Xn(p). 
isomorphism. 


(115) 

Then  <p  maps 
Hence  ip  is  an 


Q.E.D. 

If  Rjt  = {z  € V ®pW|  rank(z)  = k),  then  <^(R*)  is  the  set  of  matrices  of  rank 
k in  Fm(p)xn(p).  In  view  of  the  isomorphism,  any  linear  map  g of  V ®pW  into  itself 
can  be  considered  as  a linear  map  of  Fm(p)x„(p)  into  itself. 

Since  Homn(p)(V,  W)  ps  Fm(p)x„(p)  by  Corollary  1 of  Theorem  2.2.2,  combining 
with  the  above  theorem,  we  derive  the  following. 


Theorem  3.3.4.  Hom„(p)(V,W)  « V ®pW. 


CHAPTER  4 

THE  p-PRODUCT  FAST  FOURIER  TRANSFORM 

In  this  chapter,  we  use  the  p-product  to  develop  a novel  formulation  of  the  Fast 
Fourier  transform,  which  we  shall  refer  to  as  the  p-product  Fast  Fourier  Transform 
(pFFT).  The  efficient  implementation  of  pFFT  employs  the  tensor  product  formulation 
of  the  Fourier  transform  introduced  in  [27],  as  well  as  the  p-product  identities  proved 
in  Chapter  2.  We  omit  the  stride  permutation  matrix  and  the  tedious  effort  associated 
with  readdressing.  This  facilitates  high  performance  on  super-computers  where  the  I/O 
cost  predominates. 

4.1.  Introduction 

The  tensor  product  offers  a natural  language  for  representing  digital  signal  process- 
ing (DSP)  algorithms.  Closely  associated  with  tensor  products  are  a class  of  permutations 
called  stride  permutations.  Since  we  induce  the  p-product  Fourier  transform  from  its 
tensor  product  expansion,  we  first  define  the  stride  permutation  matrix  and  several  basic 
identities  between  the  tensor  product  and  the  p-product.  Such  identities  will  be  employed 
as  theoretical  building  blocks  in  subsequent  sections. 

Definition  4.1.1.  Suppose  that  N = rs.  The  N-point  stride  s permutation  matrix  P( N,  s) 
is  defined  by 

P(/V,  s)(x  <8>  y)  = y <8>  x,  (116) 

where  x and  y are  arbitrary  vectors  of  sizes  r and  s,  respectively. 
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This  is  to  say  that  the  action  of  P (N,s)  on  an  arbitrary  vector  x of  size  N is 


P(N,s)x.  = P(N,s)[x0,xi,...,xn-i]' 

[xq,  X2 si  • . • ^(r— l)s5  x\i  ®s+l j • • • i *^(r— l)s+l  j ®2»  • • • 5 xs—li  • • • j xrs— l]  • 


(117) 


Definition  4.1.2.  For  any  a 6 F/Xm,  we  define  co/(a),  the  column  vector  of  a,  as 

co/(a)  = [all5  a12, . . . aim,  a2i,  a22,  • • • , a2m,  • • • , a/i,  a/2, . . . , a/m]',  (118) 

and  the  row  vector  of  a as  row( a)  = (co/(a)/. 

We  now  state  and  prove  some  basic  identities  concerning  the  /^-product,  tensor 
product,  permutation  matrices,  and  Fourier  matrices. 

Lemma  4.1.3.  Suppose  N = rs  and  a e Frxs,  then 

col{s! ) = P(iV,s)co/(a).  (119) 


Proof:  Let  x = co/(a)  = [a:0,  x\, . . . , zrs_1]\  then 


a = 


3?0  X\  ...  Xg — J 

XS  *£$4-1  • • • *^25—1 


(120) 


-x(r— l)a  ^(r— l)s+l  •••  xrs— 1. 

By  Definition  4.1.1,  we  have 
P(7V,  s)x  = P(7V,  a)[x0,  xi,  • • • , ®jv-if 

[^0?  xs>  x2si  • • • x(r—l)si  xl>  *^3+1 1 • ■ • ? x(r—  l)s+l ix2 1 • • • i xs— 1 , 

= [yo,yi,...,ys-i]\ 

where 

y k — xk+si  • • • , xk+(r—  l)s]  • 


(121) 


(122) 


Moreover, 
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Therefore, 


II 

Xo 

X\ 

x 3 

XS+1 

• x(t  — 1)s 

x(r— l)s+l 

7 

H 

x2s— 1 

1 

7 

fc. 

Co/( a ) [^0)  Xs,  • • • x(r— l)s5  • • • i xs—  1)  • • • ®ra— l] 

= [yo,yi,..  • ,ys-if  = P(iV,s)x  = P(iV»co/(a). 


Lemma  4.1.4.  Let  a £ Fsxs,  with 


and  b £ Frxr,  with 


a o a\ 

as  as+1 

- a(s—  l)s  a(s— l)s+l 


bo  b\ 

br  ^r+l 

-^(r-l)r  ^(r-l)r+l 


Let  N=rs  and  x be  a vector  of  length  N. 


<*s-l 

U2s-1 

as2 


br- 1 

&2r-l 

bT  2 


(a)  //yi  = P(iV,  s)(b  @rx),  and  y2  = ro^(b)  ©rx,  then 


yi  = co/(y2). 


(b)  //y£  Faxr,  then 


(123) 


(124) 


Q.E.D. 


(125) 


(126) 


(127) 


row( a)  ©sy  = (a@sco/(y))'. 


(128) 
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Proof:  (a)  Associate  anrxs  matrix  a with  x as  follows: 

X\  . . . Xg—\ 

^5  + 1 • • • *^25  — 1 

x(r— 1)5+1  • • • xrs— 1 _ 

Thus,  by  Definition  4.1.2,  x = co/(a). 

Using  Theorem  2.2.12  (a)  of  Chapter  2,  we  obtain 


x0 

xs 

- x(r— l)s 


(129) 


b ©rx  = b ©rco/(a)  = co/(b  x a).  (130) 

Alternatively,  by  Theorem  2.2.12  (c), 

y2  = rou>(b)  @rx  = (b  x a)'.  (131) 


Combining  this  with  Lemma  4.1.3,  (a)  follows. 

Similarly,  by  using  Theorem  2.2.12  (a)  and  (b),  we  can  prove  (b). 

Q.E.D. 

More  generally,  if  x is  an  N x t matrix,  then  the  following  lemma  can  be  proved. 


Lemma  4.1.5.  Let  N=rs,  b be  the  matrix  defined,  in  Lemma  4.1.4,  and  x be  an  N x t 
matrix.  Suppose 


Yi  = P(N,  s)(b  ©rx)  and 
y2  = row(b)  0rx, 


(132) 


then  co/(yi)  = co/(y2). 

Proof:  Write  the  r x r matrix  b 


fbil 


br  J 


as 


, where  bj  [^(i— l)r>  ^(i— l)r+li  •••>  ^ir— l]- 


(133) 
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Then 


yi  = P(7V,5)(b©rx)  = P(iV,5) 


/ 

bi  ®rx' 

\ 

©rX 

V 

. br  ®rx. 

/ 

= P (N,s) 


gi 

g2 

Lgr  J 


, with  g,  £ Fsxt. 


On  the  other  hand. 


y2  = row( b)  ©rx  = [bi,  b2,  . . . , br]  @rx 


= [bi©rx,  b2©rx,  br©rx] 

= [gl,  g2,  - gr]- 

The  result  now  follows  from  Lemma  4.1.3. 


Lemma  4.1.6.  Let  a £ Fm2xm2  and  f £ Fm2.  If  a can  be  separated  as  a 
then 

g = af  = col[row{&T)  ©mco/(rou;(ac)  ©mf)). 

Proof:  By  the  identities  of  the  tensor  product,  we  have 

3.  — 3C  ® 3r  = (Im  ® ® Im ) 

= P(m2,m)(ar  <g>  Im)P(m2,  m)(ac  ® Im). 
According  to  Lemma  4.1.4, 

P(m2,m)(ac  ®Im)f  = co/(rou;(ac)  ©mf). 

The  result  follows  from  Lemma  4.1.4. 


(134) 

(135) 
Q.E.D. 

a c ® 

(136) 

(137) 

(138) 
Q.E.D. 


The  next  lemma,  which  is  a key  to  subsequent  work,  is  proved  in  [27]. 
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Lemma  4.1,7.  Let  N=rs  with  r and  s being  any  integers.  The  N-point  Fourier  transform 
can  be  factored  as 

F (N)  = (F(s)  <g>  Ir)Tr(N)P(N,  s)(F(r)  <g>  I5),  (139) 

where  P (N,  s)  is  an  N-point  stride  s permutation  matrix  and  T r(N)  is  a diagonal  matrix  as 
T r(N)  = diag  (l, . . . , 1;  1,  w, . . . , it/-1; . . . ; 1,  ws~x, . . . , (140) 

with  w — e2*l/N . 


4.2.  Basic  Properties  of  Fourier  Transform  Matrix 


The  Fourier  transform  matrix  of  order  N,  denoted  by  F(N)  is  defined  as 


F (N)=  , wN  = e^'N,i=V =T. 


The  conjugate  of  w jv,  denoted  by  w*N  is 


(wN)  = e 


* = ,(-2  *i/N)  = w-i 


N ’ 


(141) 


(142) 


and 


(^Ar)  ~ wn 

Direct  computation  shows  that 


—kmodN  _ N-k 
LLJ  \r  — W N 


- fry 


(143) 


F(N)F(N)*  = NIn.  (144) 

The  inverse  Fourier  transform  matrix  is 

F(JV)-1  = Lf(JV)',  (145) 

and  F(N)  is  symmetric,  i.e., 


F(N)'  = F(N). 


(146) 
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4.3.  An  Example  of  the  p-product  Fast  Fourier  Transform  Algorithm 
The  8-point  Fourier  transform  of  a vector  x is  given  by  the  formula 

7 

yk  = ^ wjkXj,  0 < k < 8,  w = e 2’r‘/8. 

j= o 


(147) 


In  [27],  the  tensor  product  formula  of  8-point  Fourier  transform  has  been  proved  as 

F(8)  = (F(4)  ® I2)T(I4  ® F(2))P(8, 4) 

= (F(4)  ® I2)TP(8, 4)(F(2)  <H>  I4), 


where 


T = diag(l,  1, 1,  w,  1,  w2 , 1,  u>3), 


and  P(8, 4)  is  the  8-point  stride  4 permutation. 
By  Lemma  4.1.7, 


(148) 


(149) 


F(4)  = (F(2)  (8)  I2)T4P(4, 2)(F(2)  <8>  I2). 


(150) 


Using  this  and  the  tensor  product  identity. 


(BC)®I  = (B®I)(C®I), 


(151) 


the  operation  F(4)  ® I2  can  be  factored  as 


F(4)  ® I2  = (F(2)  (8)  I4)(T4  (8)  I2)(P(4, 2)  ® I2)(F(2)  ® I4),  (152) 

where 

T4  = diag(l,  1, 1,  w2).  (153) 

Placing  the  preceding  result  into  Equation  (148),  we  have  the  factorization  of  F(8), 


(F(2)  ® I4)(T4  (8)  I2)(P(4, 2)  (8)  I2)(F(2)  ® I4)TP(8, 4)(F(2)  ® I4). 


(154) 


We  organize  the  computation  into  stages  by  setting 

Y3  = TP(8,4)(F(2)  <8>  I4), 


44 


Y2  = (T4  (8)  I2)(P(4, 2)  (g)  I2)(F(2)  ® I4),  (155) 

Yi  = F(2)  ® I4. 

Applying  Lemmas  4.1.4  and  4.1.5  as  well  as  Theorem  2.2.8  to  these  stages,  we  induce 
the  /^-product  formulae  of  the  preceding  equations. 

First,  by  Theorem  2.2.8,  for  any  8-point  vector  x,  we  have 

Y3x  = TP(8,4)(F(2)  ® I4)x 


= TP(8,4)(F(2)  ©2x). 


(156) 


According  to  Lemma  4.1.4  (a), 


p(8, 4)(F(2)  ©2x)  = co/(rou>(F(2))  02x).  (157) 

Setting  w = rou;(F(2))  = [1,1,— 1,1],  we  obtain 


p(8, 4)(F(2)  02x)  = col( w 02x).  (158) 


Let  v = P(8, 4)(F(2)  02x),  and  u = w 02x.  By  the  definition  of  the  GMP,  we 


know  that  v is  a column  vector  of  length  8 and  u is  a matrix  of  size  4x2.  If  we 


set  v = [t?o,  ui, . . . , ^7]  , then  since  v = co/(u),  we  have 


Vo 

Vi 

V2 

V3 

V4 

V5 

V6 

V7 

Hence,  we  derive  the  following: 

TP(8,4)(F(2)  ©2x)  = Tv 


(159) 


[vo 


,vi,v2,wv3,v4 


,w2v5,v6,w3 


(160) 


( 

'1 

1 ' 

’ VO 

V\ 

\ 

( 

’^0 

Vl 

\ 

1 

1 

W 

2 

* 

V2 

V3 

= col 

V2 

WV3 

2 

W 

v4 

V5 

V4 

W V5 

V 

1 

W3  _ 

_V6 

v7_ 

) 

V 

,v6 

W3V7 

J 

= col 
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Let 


It  follows  that 


D2 


'1  1 ' 
1 w 
1 w2 
1 w3 


Y3x  = TP(8, 4)(F(2)  0 I4)x  = co/(D2  * (w  02x))  = co/(y3). 


Second,  by  the  tensor  product  identity, 

(BC)0I  = (B®I)(C®I) 


we  can  write  Y2co/(y3)  as 

Y2co/(y3)  = (T4  ® I2)(P(4, 2)  ® I2)(F(2)  ® I4)co/(y3) 

= [(T4P(4, 2)(F(2)  ® I2))  0 I2]co/(y3). 
Applying  Theorem  2.2.8,  the  above  equation  becomes 

Y2co/(y3)  = [(T4P(4, 2)(F(2)  0 I2))  0 I2]co/(y3) 

= (T4P(4, 2)(F(2)  0 I2))  04co/(y3). 
Now,  using  Theorem  2.2.12  (a),  we  obtain 

Y2co/(y3)  = (T4P(4, 2)(F(2)  0 I2))  ©4co/(y3) 

= co/(T4P(4, 2)(F(2)  0l2)y3) 

= Co/(T4P(4,2)(F(2)e2y3))- 

By  Lemma  4.1.5, 

col{ P(4, 2)(F(2)  02y3))  = col( w 02y3). 


If  we  set 


VO 

Vi 

V2 

V3 

v4 

V5 

V6 

v7 

(161) 


(162) 


(163) 


(164) 


(165) 


(166) 


(167) 


v=P(4,2)(F(2)@2y3)  = 


(168) 
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then,  similar  to  the  derivation  of  the  first  stage,  we  have 

u = w ©2y3  = 


and 


Vo  V\  V2  V3 

w4  u5  v6  v7 


col(  T4v)  = col 


where 


/ 

'1 

0 

0 

0 ' 

"w0 

Vi 

\ 

0 

1 

0 

0 

Y 

v2 

V3 

0 

0 

1 

0 

A 

V4 

V5 

V 

0 

0 

0 

w2  _ 

,U6 

v7_ 

) 

/ 

■ 1 ■ 

W0 

V\  ' 

\ 

( 

wo 

V\ 

\ 

col 

1 

* 

v2 

W3 

= col 

u2 

W3 

1 

V4 

v5 

W4 

W5 

w2  _ 

.*>6 

V7  _ 

) 

\ 

W2Vo 

w2v7 

/ 

= col 


1 1 1 
1 1 w2 


D4  = 


w 


* U ^ = col{ Di  * (w  ©2y3)), 


1111 
11  w2  w2 


(169) 


(170) 


(171) 


Hence 


co/(Y2co/(y3))  = col( Di  * (w  02y3))  = co/(y2).  (172) 


Finally,  by  similar  derivation,  we  induce  the  following: 

Yico/(y2)  = (F(2)  0l4)co/(y2)  = (F(2)  ® I2)  ©4co/(y2) 


(173) 

= co/((F(2)  <g>  I2)y2)  = co/(w  02y2)  = coZ(yi). 

Since  y3  is  a row  vector,  y'j  = co/(yi). 

It  now  follows  that  for  any  8-point  vector  x,  its  p -product  Fourier  transform  formula 
is  given  by 


y = ^ff(8)x  = ^(w  ©2°i  * (w  ©2D2  * (w  ©2x)))'. 


(174) 
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4.4.  The  ^-product  Fast  Fourier  Transform  for  N = 2” 


The  A-point  Fourier  transform  of  a vector  x is  given  by  the  formula 
N- 1 

Vk=Yj  wJkxn  0 < k < N,  w = e2*,/N.  (175) 

j= o 


Theorem  4,4.1,  Let  TV  = 2".  The  tensor  product  formula  of  N -point  Fourier  transform 
is  given  by 


= ^n((F(2)<8)I2n-l)(Tt  + 1 ®I2n—l)(PI  + 1 ®I2„-i-l))^(F(2)®I2n-l), 


(176) 


where  T,  is  a 2'  x 2'  diagonal  matrix  given  by 


?;  = diag(l,l,l,w2  , l,*/;2'2"  . . l,to(2’  1)2" 


(177) 


and  P,  is  the  2,-point  stride  2'  1 permutation  matrix. 

Proof:  We  will  prove  this  by  induction.  In  the  last  section,  we  have  proved  that  this 
theorem  holds  for  n < 3.  Assume  for  M = 2n~1  we  have 

/ w— 2 \ 


F(M)  = n ((F(2)0I2.-O(T<+1<8)I2-  -•-2)(P»+i  ® I2"--2)) 


.t=i 


(178) 


(F(2)0l2»-a). 

We  prove  that  this  assumption  implies  formula  (176)  for  N = 2n. 
Since  JV  = 2M,  by  Lemma  4.1.7  we  have  that 


F(N)  = (F (M)  0 12)TM(A)P(iV,  2)(F(2)  0 Iat)-  (179) 

Replacing  F(Af)  by  Equation  (178),  the  above  equation  becomes 

W)=  (n  «F(2)  ® I2.-.  )(T.+1  ® I2n-,-,)(Pl+1  0 )))  (180) 

(F(2)  ® I2.-,)T„(N)P(N, 2)(F(2)  ® I2.-. ). 
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Since  Tm(N)  = T2.-i(iV)  ® I2o,  it  follows  that 

(n— 2 

n ((F(2)®I2»-i)(T<+10I2»  -i-l)(Pi+l  <8>I2n-i-l)) 

1=1 

(F(2)  ® I2n-1  )TJV/ (TV)P(iV,  2)(F(2)  ® I2„-,)  (181) 

n — 1 \ 

JJ  ((F(2)  ® I2n-i  )(TI+1  ® Ia-i-i  )(P,+1  ® ))  j (F(2)  ® I2-i ). 


Q.E.D. 

We  next  employ  the  preceding  theorem  to  derive  the  p -product  Fourier  transform. 
The  notation  ones («)  will  denote  a row  vector  of  size  n all  of  whose  elements  are  equal 
to  1. 


Theorem  4.4.2.  Let  N = 2".  The  p-product  Fourier  transform  of  any  N -point  vector 
x is  given  by 

y = -^F(JV)x  = ^(w  ©2Di  * (w  ©2D2  * (•  • • (D„_1  * (w  ©2x)))))',  (182) 

where  w = roiu(F(2))  = [1, 1,  -1, 1],  and  D*  is  a 2k~ 1 x 2n~k+1  matrix  associated  by 
the  diagonal  elements  of  T*,  which  is  defined  as 


Dfc  = [djbdfc]  ®ones^2n 


(183) 


with 


d * 


^ l,w,w2, 


2 iri/N 


and  w — e 


(184) 


Proof:  We  partition  the  computation  of  Equation  (176)  into  stages  by  setting 

Yb  = (T„)(Pb)(F(2)®I2-i), 
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Yn_i  = (Tn_i  <S>  l2)(Pn-i  ® l2)(F(2)  <g>  I2— 1 ), 

Y,  = (T,'  ® l2n-i)(P , ® l2n-i)(F(2)  <8>  I2n-1  ),  (185) 


Y2  — (T2  ® l2n_2 )(P 2 ® l2"-2)(F(2)  <g)  Ljn-i), 

Yl  = F(2)  <g)  I2n-1. 

In  order  to  complete  the  proof,  we  first  verify  the  following  formula. 


co/(TfcP*(F(2)  ® I2*-i)x)  = col(Dk  * (w  ©2x)),  (186) 


where  D*  is  a matrix  having  the  same  size  as  w @2x  and  is  associated  with  the 
diagonal  matrix  T*,  w = row(F(2))  = [1, 1, 1,  -1],  and  x is  a 2k  x 2n~k  matrix. 

In  fact,  if  x is  a 2k  x 2n~k  matrix,  then  by  Theorem  2.2.8  and  Lemma  4.1.5  we  have 
co/(Pt(F(2)  ® I2*-i)x)  = co/(Pfc(F(2)  ©2x)) 


(187) 


= col(row( F(2))  ©2X)  = co/(w  @2x). 

Now  let  v = P*(F(2)  <g>  I2*-i  )x,  and  u = w @2x.  By  the  definition  of  the  GMP,  we 
know  that  v is  a matrix  of  size  2k  x 2n~k  and  u is  a matrix  of  size  2k~ 1 x 2n-fc+1. 
If  we  write  v as 


o 

i 

V^n—k  — l 

Vo 

. V^n  — fc+1 J 

= 

Vl 

-v(2*-l)2"-*  • 

1 

7 

• c 
cs 
£ 

1 

rH 

1 

• Ji 

CS 

> 

1 

(188) 


where 


vt  — [u,'2»-*>l’i2n-A:+l?  • • • Ju(i+l)2n-fc-l]j 


(189) 


then,  since  co/(v)  = co/(u),  it  follows  that 
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u = w ©2x  = 


00 

02"-*  — 1 

02 

. V<2n— fc  + 1 j 

02"-*+ 1 

03-2n_1  — 1 

03-2"-* 

04. 2"-*  — 1 

-v(2k-2)2n~k 
V0  Vi 

• ^(2*  — l)2n-t  — 1 

v(2k-l)2n~k  • 

02" -1 

V2  V3 

. 

-v2*-2  v2*-l  - 

Hence, 


co/(Tfcv)  = col 


/ 

■ 1 ■ 

v0 

\ 

/ 

Vo 

101 

Vl 

101  Vl 

1 

V2 

V2 

W2 

* 

V3 

= col 

102V3 

V 

1 

V2*_2 

) 

v2fc— 2 

. w2k  — l - 

-v2*-l- 

. 0^2*  — 1 v2*  — 1 - 

= col( ( [dj, dfc]  ® ones ^2"  * ) ) * u)  = co/(D*  * (w  ®2x)), 


where 


D*  = [dJ,djt]  <8>  ones ^2"  fc), 


with 

djt  = ^[1, 101,102, . . .W2k-i_if  ^ , and  wj  — w J = e2ir,^N . 


Now  we  apply  Equation  (186),  tensor  identity  (163),  Theorem  2.2.8  and 
2.2.12  to  complete  the  proof.  If  x is  an  N-vector,  then 

Ynx  = TnPn(F(2)  (8)  l2«-i)x 


(190) 


(191) 


(192) 


(193) 


Theorem 


= col( D„  * (w  ®2x))  = co/(yB), 


(194) 


and 
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Yn-icol(yn)  = (Tn_i  0 I2)(P„_i  0 I2)(F(2)  <g)  I2n-i)co/(yn) 

= ((Tn_iP„_i(F(2)  0 I2n-2))  0 l2)co/(y„) 

= ((T„-iP„-i(F(2)  0l2n-2))  02„_ico/(yn))  (195) 

= co/(Tn_iP„_1(F(2)  0 I2n-2)yn) 

= col( Dn_i  * (w  02y„))  = co/(yn_i). 

Similarly,  for  any  integer  n — 2 > i > 1,  we  have  that 

Y,' co/ (y,_)_i ) = (T,-  0 I2r.-.)(Pj  0 I2n-i)(F(2)  0 I2»-i  )co/(y,) 

= ((T,P «(F(2)  0 12.-i ))  0 12n-i)co/(y,) 

= (T,P,(F(2)  0 12— i))  ©2.co/(yi)  (196) 

= co/(TlP,(F(2)0l2.-1)y,) 

= co/(D,  * (w  02y,))  = co/(y,). 

Thus, 

y=iF(Jf)x=iY,Y!...Y,x 

= -^co/(w  ®2(D2  . (w  ®2(. . . (D„  » (w  ®2x)))))).  (W7) 


Since  y is  a vector,  the  result  follows. 


Q.E.D. 


4.5.  The  p-product  Fourier  Transform  for  N = an. 


The  Appoint  Fourier  transform  of  a vector  x of  length  N is  given  by  the  formula 

N-l 

yk=Yl  w3hxv  0 < k < N,  w = e2*l/N.  (198) 

j=0 


By  modifying  the  proof  of  Theorem  4.1.4,  we  obtain  the  following  theorem. 
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Theorem  4.5.1.  Let  N = an  with  a > 2.  The  tensor  product  formula  of  any  N-point 
Fourier  transform  is  given  by 


F(N)  = ^n((F(a)  ® Ia— i)(T,+i  ® Ia.-i-i)(P,+i  (8)  (F(a)  ® Ia-i),  (199) 


where  Ti  is  an  a1  x a*  diagonal  matrix  given  by 


T i 


n— « 

diag(l, . . . 1;  1, w 

u,(“— 1)“"~'-  • 

. . ,tt)(“-1)(a'-1-1)a 

), 


(200) 


and  Pj  is  the  a' -point  stride  a'  1 permutation  matrix. 

Now,  using  the  above  theorem  and  substituting  2 for  a in  the  proof  of  Theorem 
4.1.5,  we  can  induce  the  p-product  Fourier  transform  of  an  N-point  vector  x. 


Theorem  4.5.2.  Let  N = a”.  The  p-product  Fourier  transform  of  any  N-point  vector 
x is  given  by 

y = jF(N)x  = -^(w  ©aDi  * (w  ©aD2  * (. . . (D„_i  * (w  ©ax)))))\  (201) 

where  w is  the  row  vector  corresponding  to  F(a), 

D,  = [d-,d-,d d“_1]  ®ones(an"'),  for  i = 1,2 - 1, 


and  d,  = ( 


1,  w,  w , . . . , w 


)' 


(202) 


is  an  a1  1 x an  1+1  matrix  associated  by  the  diagonal  elements  o/T,  with  w = e2x,/Jv. 


4.6.  The  p-product  Fourier  Transform  for  Prime  Factorization  of  N 


In  the  previous  sections  we  formulated  p-product  versions  of  the  FFT  for  signals 
of  length  N where  N = an.  However,  in  many  real  applications  N is  not  equal  to  a 
power  of  some  integer  a.  In  such  cases  it  is  possible  to  use  the  prime  factorization  of  N 
in  order  to  obtain  sparse  factorizations  of  F(iV). 
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Theorem  4.6.1.  If  N = p\p2  • • - Pn>  with  pi  being  a prime  number,  then  the  tensor  product 
formula  of  any  N -point  Fourier  transform  is  given  by 

F(i V)  = ( fl  ((F(p,)®I».)(T,+i  ®Ia.„)(P,+i  ®Ia,+,)) 'W„)®IsJ, 

Vi=i  / (203) 

N 

where  8{  — N/p{ , A,-  = , 

PlP2  • • • P. 

where  T,+i  is  a p\p2  . . . pI+i  x p\p2  . ..pi+ 1 diagonal  matrix  as  follows: 

T«'+i  = Tp.+i  (P1P2  • • • Pi+l) 

= (diag  (1, . . . , 1;  1,  tu, . . . , wPi+1_1; . . . ; (204) 

l,wpip2-p'~1, . . . , w(i ’-h-IXpi-p;-!)))^ 5 
and  P,-+1  is  a p\p2  . . -Pi+i— point  stride  p\p2  . ■ - Pi  permutation  matrix. 

Proof:  We  prove  Theorem  4.6.1  by  induction.  By  Lemma  4.1.7,  we  know  Equation 
(203)  is  valid  for  n= 2.  Assuming  that  I<  = • • -Pn-i.  Equation  (203)  holds,  we 

prove  this  for  N — p\p2  . . .pn-  Since  N = p\p2  . . .pn  = Kpn,  apply  Lemma  4.1.7 
once  more  in  order  to  obtain 


F(iV)  = (F (K)  0 IPn)TPn(N)P(N,pn)(F(pn)  0 1^).  (205) 


By  induction  hypothesis,  we  have  that 

(n— 2 

n (F(k)  ® I».)(Ti+l  ® Ij,t,)(P.-+l  ® Ia,„)  ) (F(Pn-i)  ® 


(206) 


with  T,,  P,-,  Ati  and  8t  defined  by  Equations  (204)  and  (203),  respectively.  Combining 
two  Equations  (205)  and  (206),  the  result  follows. 


Q.E.D. 

Using  the  preceding  theorem  and  the  p-product  identities,  modification  of  the  proof 
of  Theorem  4.1.5  yields  the  following  theorem. 
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Theorem  4.6.2.  Let  N = p\p2  . . .pn,  with,  each  pt  a prime  number  and  pi  < pi+i-  The 
p-product  Fourier  transform  of  any  N -point  vector  x is  given  by 

y = = ^(W1  ©PlDi  * (w2  ©P2D2  * (•  • • (D»_1  * (w„  ©p„x)))))',  (207) 

where  w,  = rou;(F(pj)), 

d°,  d,- , . . . , d^,_1  ® ones(A,), 

\f' 

is  a pip2  . . .pi- 1 x piPi+ 1 . . .pn  matrix  associated  by  the  diagonal  elements  ofTi  and 
w = 


Di  = 


with  d 


< = ([>.■ 


*2  — 
, . . . , 


Y)  ’A*  = N/PlP2  ■ • • Pi  i Ti  =PlP2  ■■■Pi-1, 


(208) 


4.7.  Two-dimensional  p-product  Fourier  Transform 

Since  the  Fourier  transform  is  separable,  it  is  easy  to  induce  the  two  dimensional 
p-product  formulation,  given  the  previous  theorems.  In  order  to  do  so,  we  need  to  define 
a function  ipn  : Fm  — i ► F/Xn  by 


a\ 

a2 

■ an 

V’n([ai,a2,...,am])  = 

an+l 

an+2 

• a2n 

, (209) 

-a(l-l)n+l 

a{l-l)n+2  ■ 

• aln  . 

with  m-ln,  and  then  apply  this  function  after  performing  row  and  column  transforms. 


Theorem  4,7.1.  Let  N — p\p2  . . . p„  and  M — q\q2  ...  qm.  For  a two  dimensional  NxM 
array  x,  its  Fourier  transform,  denoted  by  y,  is  given  by 

y = ~F(N)xF(M),  (210) 

and  its  corresponding  p-product  factorization  is  given  by 

y = W!  ©PlDi  * (w2  ©P2D2  * (. . . (Dn_j  * (w„  ©Pnx'))))),  (211) 

with 

x = F (M)x'  = il>N[v i ©giEi  * (v2  ©g2E2  * (. ..  (Em_i  * (vm  ©9mx'))))],  (212) 
and  w,,D,  and  v,-,Ei  are  defined  in  Theorem  4.6.2. 
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4.8.  Computational  Cost 


The  number  of  arithmetic  operations  required  to  carry  out  a computation  is  an 
important  part  of  the  cost  of  the  computation  and  has  traditionally  attracted  much  attention. 
The  computations  required  in  the  p-product  Fourier  transform  are  of  the  same  order  as 
the  Cooley-Tukey  FFT;  i.e.  0(7V  log  N)  operations.  However,  the  reordering  process 
required  in  the  customary  FFT  imposes  a performance  penalty  upon  modem  architecture 
machines.  H.  S.  Stone  [25]  has  shown  that  up  to  0(N)  time  is  required  to  execute  the 
reverse  operation  for  N being  a power  of  two.  Thus,  the  p-product  Fourier  transform 
algorithm  can  save  0(A)  time  by  eliminating  reordering. 

We  have  implemented  the  pFFT  on  a sequential  machine  (SUN  3 Workstation) 
using  Matlab  3.0  software.  As  shown  in  Table  1,  the  time  required  for  the  new  algorithm 
is  less  than  the  time  required  by  the  Cooley-Tukey  method. 

Although  the  pFFT  exhibits  similarities  to  the  Stockham  autosort  algorithm  [27], 
it  has  the  advantage  of:  a)  concise  expression,  by  replacing  the  permutation  matrix  and 
the  tensor  product  with  the  single  p-product,  and  b)  a standard  implementation,  given 
the  p-product  as  a basis  operation. 


Table  1 Time  required  to  compute  Fourier  transform  of 
A-vectors  on  the  SUN  3 Workstation  using  Matlab  3.0  with 
Cooley-Tukey  FFT  and  the  p-product  Fourier  transform  algorithms. 


N= 2 

AM 

N= 8 

AM  6 

N=32 

AM4 

AM  28 

N=256 

AM  12 

AM  024 

0.018 

0.033 

0.047 

0.063 

0.086 

0.122 

0.177 

0.279 

0.491 

0.909 

tc 

0.032 

0.048 

0.066 

0.090 

0.116 

0.167 

0.230 

0.369 

0.592 

1.110 

t c'  time  required  for  using  Cooley-Tukey  FFT  algorithm.  tp:  time  required  for 
using  p-product  Fourier  transform  algorithm. 


CHAPTER  5 

THE  ^-PRODUCT  WALSH  TRANSFORMS 


The  Walsh  transform  and  the  generalized  Walsh  transform  are  similar  to  the  Fourier 
transform  but  exhibit  simpler  formulations.  The  Walsh  transform  was  defined  in  1923 
by  J.L.  Walsh  [28],  although  Hadamard  (1893)  [14]  achieved  a similar  result  by  the 
application  of  certain  orthogonal  matrices,  generally  called  Hadamard  matrices,  which 
contain  only  the  entries  +1  and  -1.  A special  set  of  these  matrices  can  be  shown  to 
be  directly  related  to  the  Walsh  series  through  a tensor  product  operation  on  a basic 
Hadamard  matrix. 

In  1931,  R.  E.  A.  C.  Paley  provided  an  entirely  different  definition  of  the  Walsh 
transform,  which  is  the  one  used  most  frequently  by  mathematicians  [11]  and  is  the 
one  used  in  this  dissertation.  This  definition  retains  the  advantages  of  the  earlier 
formulations  (good  convergence  properties  for  2nth  partial  sums  and  close  analogy 
with  the  trigonometric  functions)  while  adding  two  new  features.  By  using  products 
of  Rademacher  functions,  Paley  provided  a vehicle  to  obtain  explicit  representations 
of  kernels  of  certain  operations  associated  with  Walsh-Fourier  series  (for  example, 
the  sequence  of  partial  sums),  and  forged  a link  between  the  theory  of  Walsh  series 
and  probability  theory.  However,  transform  algorithms  derived  from  this  formulation 
require  either  that  the  input  data  vector  is  arranged  in  bit-reversed  order  of  placing, 
or  that  a similar  rearrangement  is  carried  out  on  the  output  vector.  This  bit-reversed 
natural  order  for  the  coefficients  is  not  convenient  for  most  signal  processing  and 
communication  applications  due  to  difficulties  of  calculation  and  interpretation.  For 
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this  reason,  programmers  do  not  prefer  this  formulation  as  much  as  the  other  two.  This 
may  also  be  the  reason  that  although  mathematicians  have  established  a very  strong 
mathematical  background  for  this  newer  formulation,  practitioners  still  do  not  apply 
them  in  industrial  problems. 

We  previously  showed  that  the  p-product  can  be  used  to  compute  the  FT  without 
requiring  the  reordering  process.  Given  similarities  between  the  Fourier  functions  and 
Walsh  functions,  we  can  obtain  the  same  results  for  Walsh-Paley  functions.  Thus,  we 
retain  the  mathematical  advantages  of  Paley’s  equation,  and  achieve  the  computation  in 
proper  order  without  requiring  pre-  or  post-ordering. 

In  1955,  Chrestenson  [4]  extended  Paley’s  idea  and  formed  a class  of  generalized 
Walsh  functions  which  include  the  Walsh  functions  as  its  special  case.  Different  from 
the  Walsh  transform,  this  definition  is  based  on  finite  products  of  Rademacher  functions 
of  order  a and  is  related  to  an  a-nary,  a > 2,  decomposition  of  the  index  to  the  function 
so  that  the  order  is  referred  to  as  an  a-nary  ordered  or  natural  series.  Generalized  Walsh 
functions  have  mathematical  properties  similar  to  those  of  Walsh  functions. 

In  this  chapter  we  discuss  the  unification  of  transform  algorithms  derived  from 
discrete  Walsh  functions  and  discrete  generalized  Walsh  functions  in  terms  of  the  p- 
product  language.  The  p-product  provides  an  important  mathematical  tool  to  combine  the 
two  kinds  of  functions  under  one  common  parasol.  One  standard  algorithm  can  carry  out 
the  computations  of  both  functions.  It  also  provides,  for  the  first  time,  a fast  method  for 
computing  generalized  Walsh  functions.  Furthermore,  we  extend  Chrestenson’s  idea  and 
use  the  similarities  and  differences  between  the  Fourier  functions  and  the  Walsh  functions 
to  derive  a representation  of  the  generalized  Walsh  functions  for  the  case  of  N being  the 
factors  of  prime  numbers.  Finally,  we  generate  a class  of  p- product  matrix  transforms 
obtained  from  corresponding  tensor  matrix  product  transforms.  This  demonstrates  that 
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any  orthogonal  or  nonorthogonal  transforms  obtainable  from  a series  of  matrix  tensor 
products  can  be  expressed  concisely  in  terms  of  the  p -product,  and  may  be  efficiently 
implemented  on  digital  computers. 

5.1.  Introduction 

To  define  the  generalized  Walsh  functions,  Chrestenson  used  the  following 
Rademacher  functions  of  order  a. 

Let  a denote  a fixed  integer  with  a > 2,  and  let  u = e27r,/“. 

Definition  5.1.1.  The  Rademacher  functions  of  order  a are  defined  by 

<f>o(x)  = a>k  if  k / a < x < (k  + 1) /a,  k = 0, ...  ,a  — 1,  (213) 

and,  for  n > 0, 

4>n(x  + 1)  = (f>n{x)  = <j)o(anx).  (214) 


Thus,  for  q=2  a Rademacher  function  of  index  n is  a train  of  rectangular  pulses 
with  2n  cycles  in  the  half  open  interval  [0,1),  taking  the  values  +1  or  -1  (see  Figure 
1).  Its  period  is  1. 

Using  Rademacher  functions  we  define  the  Walsh  functions  of  order  a as  follows. 

Definition  5.1.2.  The  Walsh  functions  of  order  a are  defined  by 

i’o(x)  = 1 (215) 

and  if  n = a\ctni  + . . . + amanm,  where  0 < aj  < a , and  n\  > ri2  > . . . > nm,  then 

Mx)  = K\(x)---<l>anZ(x)- 


(216) 


59 


<t>0  (*> 


1 

o. 

-i 


01  to 


1 

0 

-1 


02  (x) 


1 


0 1 
Figure  1.  A set  of  Rademacher  functions 

For  convenience  we  let  ipa  denote  the  set  of  Walsh  functions  of  order  a.  Thus, 

2 denotes  the  set  of  functions  defined  by  Walsh.  For  a > 2,  %j)a  is  referred  to  as 

the  set  of  generalized  Walsh  functions.  It  has  been  proved  that  is  orthonormal  and 

complete  in  L(l,  0)  [4], 


5.2.  The  p-product  Walsh  Transform 
5.2.1.  Basic  Properties  of  the  Walsh  Transform 

Recall  that  discrete  Walsh  functions  are  sampled  versions  of  the  set  of  continuous 
Walsh  functions.  For  simplicity,  we  first  consider  Walsh  functions  of  order  a=2. 

Sampling  of  the  previously-defined  Walsh  functions  at  N = 2"  equidistant  points 
results  in  an  N x N matrix.  We  denote  such  matrices  by  Hw(JV),  since  they  can  be 
obtained  by  reordering  the  rows  of  Hadamard  matrices. 

The  element  h(x,u)  of  Hw(iV)  can  be  generated  by  setting 

n— 1 

h(x,u)  = JJ 
»=1 


9 


(217) 
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where  bk(z)  is  the  kth  bit  in  the  binary  representation  of  z.  For  example,  if  n = 3 and 
z = 6 (110  in  binary),  we  have  that  bo(z ) = 0,  h(z)  — 1,  and  62(2 ) = 1- 
Direct  computation  shows  that 

H w(N)  x UW(N)  = Nljf.  (218) 

Therefore,  the  inverse  Walsh  transform  matrix  is  Hiy(iV)  itself  and  H^(jV)  is  symmetric; 
i.e. 

Hw(N)  = HW(N)'.  (219) 


The  Walsh  transform  of  an  iV-point  real  array  f(x),  denoted  by  g(u),  is  given  by 

g(u)  = jf  f(x)h(x’  u)  = /(x)  II  (-l)6'W6’*-1-l(u).  (220) 

x=0  x—0  i=0 


5.2.2.  An  Example  of  the  p-product  Walsh  Transform  Algorithm 

The  continuous  Walsh  functions  and  their  corresponding  matrices  are  shown  in 
Figure  2,  for  the  case  of  N= 8.  By  rearranging  (in  bit-reverse  order)  the  rows  of  the  8 x 8 
Walsh  matrix,  we  obtain  an  8 x 8 Hadamard  matrix  shown  in  Figure  3b.  Figure  3a  is  its 
corresponding  Hadamard-ordered  continuous  Walsh  functions. 

The  Walsh  transform  can  be  computed  by  a fast  algorithm  identical  in  form  to  the 
successive-doubling  method  given  for  the  FFT.  The  only  difference  is  that  all  exponential 
terms  uff  are  set  equal  to  one  in  the  case  of  the  fast  Walsh  transform  (FWT). 
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1 

Vo  0 


Vi 


v2 


v3 


V4 


V5 


V6 


O 174  1/2  3/4  1 x*' 


a 


Hw(8)  = 


1 

1 

1 

1 

1 

1 

1 

1 


1 1 

1 1 

1 -1 

1 -1 

-1  1 

-1  1 

-1  -1 

-1  -1 


11111 
1 -1  -1  -1  -1 
-11  1-1  -1 
-1  -1-111 
-11-11  -1 
-1  -1  1-11 
11-1-11 
1-1  1 1 -V 


b 


Figure  2.  Walsh  functions  for  N=  8. 

a.  Continuous  Walsh  functions,  /V=8;  b.  Discrete  Walsh  functions,  /V=8. 
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. 1 

walh(0,x^ 

-1 


walh(1,x) 


walh(2,x) 


walh(3,x) 


walh(4,x) 


walh(5,x) 


11111111 
1-1  1-1  1-1  1-1 

1 1-1-1  1 1-1-1 


walh(6,x) 


Hw<8> 


1-1-1  1 1-1-11 
1 1 1 1 -1  -1  -1  -1 
1-1  1-1-11  -11 


walh(7,x) 


11-1-1-1-1  11 

1-1-11-11  1-1 


0 T74  T72  374  1 **  x 

a 


b 


Figure  3.  Hadamard-ordered  Walsh  functions  for  N= 8. 


a.  Hadamard-ordered  continuous  Walsh  functions,  7V=8; 


b.  Hadamard-ordered  discrete  Walsh  functions,  7V=8. 
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Thus  for  iV=8,  we  can  decompose  Hw(8)  as 


HW(8) 


+ +00000  0' 
+ -000000 
00++0000 
00  + -0000 
0000  + + 00 

0000  + -00 

000000+  + 
000000  + 


X 


+ 

0 

+ 

0 

0 

0 

0 

0 


0 

+ 

0 

+ 

0 

0 

0 

0 


+ 

0 

0 

0 

0 

0 

0 


0 

+ 

0 

0 

0 

0 

0 


0 

0 

0 

0 

+ 

0 

+ 

0 


0 

0 

0 

0 

0 

+ 

0 

+ 


0 

0 

0 

0 

+ 

0 

0 


0 

0 

0 

0 

0 

+ 

0 


'+  0 0 0 + 0 0 0' 

0 + 000  + 00 

00  + 000  + 0 

000  + 000  + 

+ 000-000 
0 + 000-00 
00  + 000-0 

. 0 0 0 + 0 0 0 -. 


(221) 


where  “+”  represents  1,  represents  -1,  and  P is  a permutation  matrix  which  performs 
the  reordering. 

As  discussed  in  the  case  of  the  Fourier  transform,  the  bit-reverse  ordering  required 
in  the  FWT  makes  the  algorithm  complicated  and  difficult  to  implement.  However,  by 
using  the  /^-product  formulation  of  the  Walsh  transform,  an  efficient  and  fast  algorithm 
can  be  derived. 

We  utilize  the  fact  that  the  p-product  Walsh  transform  expression  equals  the  FWT 
decomposition.  Rather  than  present  the  general  proof,  we  outline  the  derivation  for  the 
special  case  N= 8.  Any  one  familiar  with  the  FWT  decomposition  should  be  able  to 
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extend  the  work  presented  here  to  the  general  case,  including  input  dimensions  other 
than  powers  of  2. 


Example  5.2.2. 1.  The  Walsh  transform  of  an  8-dimensional  vector  function  f = 

[/(0), /(l),..., f(7)]'  can  be  expressed  as 

g = ^Hw(8)f  = i[wi  ©2(w2  ©2(w3  ®2f))]\  (222) 


where  wt  = [1,  1,  1,  —1]  for  i= 1,  2,  and  3. 

Proof:  Let  N = 23  = 8.  The  Walsh  transform  of  f is 


1 

g(^)  = tv  ^ f(™)h(n,m),  n = 0,1,. ..,iV  - 1.  (223) 

m—0 


By  modifying  the  Cooley-Tukey  method.  Equation  (223)  can  be  decomposed  as 

g(i,A,Jb)  4t  (-1)**“  E (-1 *)'111  <224> 


ko= 0 


*1=0 


*2=0 


where 

n = 4j2  + 2ji  + j0,  h, it, jo  = 0 or  1, 

m = 4&2  + 2&i  + &o,  ^2,  ^i,  ^o  = 0 or  1. 
Define  the  intermediate  arrays 


(225) 


Ai(jo,  ko)  = (-l)**2/(fc2,*i,fco),  (226) 

*2—0 


1 

A2(jo, ji,  &o)  = X]  (-lf^AjOo,*!,^),  (227) 

*i=0 

and 

l 

As(jo , ji , j2)  = (-l)j2A:oA2(jo,  ji, &o)- 

*o=0 


(228) 


65 


Let  a = W3  ®2f,  b = w2  ©2a,  and  c = wj  ©2b.  By  definition  of  the  p-product, 


aji  = 22  w3(ip  + k2)f(k2N(p)  + j),  for  0 < j < N(p),  0 < i < 2, 
£2=0 

1 N(p) 


ba(0,t)  = J2  Wl(PP+kl)a(klZp+a)V  for  0 < a < 2 » 

0</3<2,  0 <t  < 2, 


*1  =0 


£(„,„)  = 22  wl(uP  + k0)bkoV,  for  0 < u < 2,  0 < v < 4, 

*o=0 

and  v = 2/3  + t.  If  we  write  cn  = C(uvy  then 


(229) 


n = 4u  + 2/?  -f  t,  with  u , /3,t=  0 or  1.  (230) 

Using  the  new  index  n,  cn  becomes  c^u  p ty  Combining  the  above  equations,  we  obtain 

^ n — C(u,v)  — C(u,0,t) 


= w^up  + ko)  *22  w2(Pp  + k 1)  22  w3 (fP  + k2)f(k2N(p)  + 

*o=0  ki=0  k-i—O  ' 


(231) 


Since  N= 8,  and  p=2,  N(p)  = y = 4.  Thus, 


Ak2  + 2ki  + &o  = m,  with  k2,k\,ko  = 0,  or,  1.  (232) 

Let  Wi(up+  ki ) = ( — l)Mil.  Equation  (231)  becomes 
1 1 1 

C(„A1)  = E (-1)"40  E (-1)"4'  E (-1)“i/(4i2  + 24,  + to),  (233) 

*o=0  *i=0  *2=0 

which  is  identical  to  Equation  (224)  with  the  exception  of  the  constant  term  where 

u = j2,(i  = ji,t  = jo.  However,  the  intermediate  arrays  are  defined  differently.  Using 

the  definition  of  the  p -product,  we  define 

1 

a(kuko,j0)  = J2  (— 1)J°*2/( 4A;2  + 2 h+h), 

*2=0 

1 

b{h>,ji,jo)  = (“ ■*■)  G(*i,*o,jo)’ 

*i=0 
1 

C{h,h,Jo)  = ^ 2 ^ b(h>,ji,jo)’ 

*o=0 


(234) 
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while  the  final  result  c^t]uj0)  contains  all  the  coefficients  of  the  8-length  Walsh 
transform  and  has  the  same  order  as  g(j2,  ji,  jo)  = g(rc).  Hence 

S = ^FHw(8)f  = jfc'  = ©2(w2  ©2(w3  ©2f))]'-  (235) 


Q.E.D. 

Observe  that  final  order  of  the  FWT  is  in  bit-reversed  form  but  the  p-product  Walsh 
transform  is  not,  which  makes  the  algorithm  more  efficient.  On  the  other  hand,  since 
wi  = W2  = W3,  all  three  stages  are  identical.  As  a result,  the  same  stage  can  be  used 
sequentially  three  times  in  order  to  carry  out  the  whole  transformation.  This  makes  the 
algorithm  more  concise  and  faster  than  the  regular  FWT. 

Figure  4 shows  the  signal  flow  graph  for  the  8-length  p-product  Walsh  transform. 
Its  similarity  to  the  FWT  is  obvious. 


/ 


a 

•0.1) 


b 

b(l.l) 


c =W 
cO.l) 


c(1.5) 

c0.3) 

cO,7) 

c(U) 

cO,6) 

c(1.4) 

c(l,8) 


+ 1 


T 


Figure  4.  Signal  flow  graph  of  8-length  discrete  p-product  Walsh  transform  of  f. 
Multiplying  factor  are  +1  and  -1,  as  indicated  by  solid  and  dotted  branches,  respectively. 
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5.2.3.  The  ^-product  Walsh  Transform  for  N = 2m 


Notice  that  the  difference  between  Equation  (174) 


y = = ®2Dl  * ®2D2  * (W  ©2X)))'> 


and  Equation  (235) 


g = -^Hw(8)f  = jc'  = ^[wj  ©2(w2  @2(W3  ©2f))]' 


is  the  multiplication  of  the  matrix  D,  after  each  stage.  In  fact,  it  has  been  proved  [29] 
that  in  the  tensor  product  factorization,  the  Fourier  transform  equation  differs  from  the 
Walsh  transform  equation  by  the  introduction  of  a scalar  (or  a matrix)  multiplication  at 
each  stage.  Such  is  the  case  for  the  p -product  expression  for  N=  8.  Since  all  the  p-product 
expansions  of  the  Fourier  and  Walsh  transforms  can  be  induced  from  the  tensor  product 
factorization,  we  can  extend  this  idea  to  the  general  case  and  induce  the  p-product  Walsh 
transform  for  N = 2m. 

First,  using  the  similarities  of  the  Fourier  transform  and  the  Walsh  transform,  we 
induce  the  tensor  product  factorization  of  the  /V-point  Walsh  transform  from  Theorem 
4.4.1. 

Theorem  5.2.3. 1.  Let  N — 2”.  The  tensor  product  formulation  of  the  N -point  Walsh 
transform  is  given  by 


where  Pt  denotes  the  2*  -point  stride  2*  1 permutation  matrix. 


Second,  similar  to  the  proof  of  Theorem  4.4.2,  we  obtain  the  following  result. 
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Theorem  5.2.3.2.  Let  N = 2m.  The  Walsh  transform  of  any  N -length  function  f can  be 
expressed,  in  terms  of  the  p-product  as  follows: 

g = ^Hw(N)f  = ^Wl  ©2(W2  02(-  • • (wm  02f)))]'.  (239> 

where  v/j  = [1,  1,  1,  —1],  for  j = 1,2,  • • • , m. 

We  call  w = [1,  1,  1,  -1]  the  core  matrix  of  the  Walsh  transform  for  input 
dimensions  which  are  powers  of  2. 

Note  that  the  p-product  Walsh  transform  Equation  (239)  deals  only  with  values  +1 
and  -1.  There  is  no  multiplication  involved.  Meanwhile,  the  reordering  process  required 
in  most  Walsh  transform  algorithms  has  been  eliminated.  Thus,  computation  of  the  new 
p-product  Walsh  transform  is  more  efficient. 

5.3.  The  p-product  Generalized  Walsh  Transform 
5.3.1.  Basic  Properties  of  The  Generalized  Walsh  Transform 

Even  though  one  can  still  use  the  Cooley-Tukey  method  to  formulate  the  FWT  for 
the  case  of  a > 2,  the  implementation  on  a computer  is  more  involved  since,  in  general, 
bit-reversal  is  a time  consuming  operation.  This  limitation  restricts  FWT  algorithm  to  a 
base-two  format  only.  However,  reordering  is  not  problematic  in  the  p-product  Walsh 
transform.  Thus,  formulation  of  the  p-product  Walsh  transform  using  integer  bases  greater 
than  two  is  as  simple  as  in  the  base-two  format.  In  what  follows,  we  show  how  to  use 
the  p-product  to  express  the  Walsh  transform  when  a is  an  integer  greater  than  2. 

When  a> 2,  sampling  of  the  continuous  generalized  Walsh  functions  at  N = an 
equidistant  points  results  in  an  N x N matrix.  We  denote  such  matrices  by  Hha(TV),  in 
order  to  be  consistent  with  the  case  of  a=2. 
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The  element  h(x,u ) of  Hw(N)  can  be  generated  as  follows: 

n — 1 

h(x,  u)  = JJ  wb'WK-'~'(u\  (240) 

i=l 

where  w = e2’n/a  and  b^z)  is  the  kth  bit  in  the  a-nary  representation  of  z.  For  example, 
if  a = 3,  N = 33,  and  z = 6 (020  is  trinary),  then  we  have  bo(z)  = 0,6i(z)  = 2,  and 
b2(z)  = 0. 

The  conjugate  of  w,  denoted  by  w*,  is 


w*  = e~2*i/a  = w~\ 


(241) 


and 


„ —kmoda 
= W 


= W 


a—k 


Direct  computation  shows  that 


H^(iV)  x HW(N)*  = NIn, 


Thus  the  inverse  generalized  Walsh  transform  matrix  is 

Hw(iV)-1  = 4Hh,(JV)*, 
and  Hw(N)  is  symmetric,  i.e., 

HW(N)' = HW(N). 


(242) 


(243) 


(244) 


(245) 


For  an  Appoint  real  array  f(x),  its  generalized  Walsh  transform,  denoted  by  g(u),  is 
given  by 

9(u)  = f(x)h(x,u)  = J2  /(*)  II  u,6<(*)Jn_<"l(*) 

x=0  x—0  i=0 


(246) 
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5.3.2.  The  p-product  Walsh  Transform  for  N = am 


Let  N = a171.  Similar  to  the  decomposition  of  the  case  a=2,  the  Appoint  generalized 
Walsh  transform  can  be  factored  as 


9(u)  = j0) 


a— 1 a— 1 a— 1 /m  — 1 \ 

= — ^ i cim~lko  ^2  wjrn~2kl  ...  ^ w3okm~1f I y;  a'kj  I , 

ko=0  fci=0  fcm_i=0  \ i=0  / 


(247) 


m— 1 

where  u — a'ji- 
i= o 

By  using  the  similarities  and  differences  between  Fourier  and  generalized  Walsh 
transforms,  and  the  derivation  methods  of  Theorems  4.5.1  and  4.5.2,  we  obtain  the 
following  two  theorems. 


Theorem  5.3.2. 1.  Let  N = am.  The  tensor  product  factorization  of  the  N -point 
generalized  Walsh  transform  is  given  by 

H w(N)  = (n  ((F(a)  <g>  Ian-i)(Pl+1  <g>  Ian-i-i))J  (F(a)  Ian-i),  (248) 


where  Pi  is  the  a1 -point  stride  a*  1 permutation  matrix 


Theorem  5.3.2.2.  Let  N — am.  For  an  N -length  function  f,  the  p-product  generalized 
Walsh  transform  can  be  expressed  as: 

g = = ^[wi  ©a(w2  ©a(.  • • (wra  ©af)))]',  (249) 

where  w i,  the  row  vector  corresponding  to  F(a),  is  called  the  core  matrix  of  the 
generalized  Walsh  transform  for  dimensions  which  are  powers  of  a. 
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5.4.  The  p-product  Walsh  Transform  for  Prime  Factorization  of  N 


It  is  convenient  to  first  derive  a formula  for  the  special  case  N=2M. 

Suppose  N=2M,  modifying  the  factorization  of  the  Fourier  transform  given  in  [7], 
the  iV-point  Walsh  transform  can  be  factored  as  follows: 

M—l  ( 1 ^ 

g(j)  = g{jiJo)  = /(n  1M  + n0)wJ2oni  \wJ^no,  (250) 

ni=0  v rii  =o  J 

where  j = 2ji  + j0  for  j0  = 0, 1 and  j\  = 0, 1, . . . , M - 1. 


Using  the  method  in  [27],  we  can  write  above  formula  in  the  tensor  product  form 


Hi v(N)  = (F(M)  0 I2)(Im  ® F(2))P(iV,  M) 
= (F(M)  0 12)P(iV,  M)(F(2)  0 IM). 


(251) 


Applying  Hn/(Ar)  to  an  A-length  function  f,  by  Theorem  2.2.8,  it  follows 


Hw(N)f  = F(M)  M)( F(2)  ©2f)).  (252) 


By  applying  Lemma  4.1.4,  we  obtain  the  p-product  Walsh  transform  of  the  function  f as 


Hw(N)i  = (row(F(M))  ®M(row( F(2))  ©2f))'.  (253) 


Note  that  M is  prime,  then  we  have  F(M)  = H w(M).  Thus,  the  preceding  equation 
becomes 

Hw{N)f  = ( row(Hw(M ))  ®M(row(Hw(2))  ©2f))'.  (254) 


Carrying  out  the  same  procedures,  and  using  Theorems  4.6.1  and  4.6.2,  the  next 
theorem  can  be  easily  proved  by  induction. 
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Theorem  5.4.1.  Let  N = p\p2  . . .pn  with  pi  a prime  number  and  pi  < p,+i.  The  tensor 
product  factorization  Hw(iV)  is  given  by 

Hiv(JV)  = ( n (F(n)®It,)(P,+i  ®lA,tl) 

V.=i  / (255) 

N 

with  6i  = N/pi,  A ,•  — , 

P1P2  ■■■Pi 

where  P,  is  p\P2  ■ ■ .pi-point  stride  pip2  • • - Pi- 1 permutation  matrix. 

Theorem  5.4.2.  Let  N = p\p2  . . . pn  with  pi  being  a prime  number  and  pi  < pl+\ . The 
p-product  generalized  Walsh  transform  of  any  N -point  vector  f is  given  by 

8 = = ^(wi  ©P1  (w2  ©P2  (. . . (wn  ©Pnf))))',  (256) 

where  w,  is  the  row  vector  corresponding  to  F(pt). 

Note  that  the  p-product  generalized  Walsh  transform  algorithm  exhibits  of 
0(N  log  N)  complexity  without  reordering.  This  is  less  burdensome  than  the  traditional 
method  which  requires  N 2 operations. 

5.5.  The  Two  Dimensional  p-product  Walsh  Transforms 

As  in  the  case  of  the  Fourier  transform,  both  the  Walsh  and  generalized  Walsh 
transforms  are  separable.  Therefore,  in  order  to  express  the  two-dimensional  Walsh 
transform  in  terms  of  the  p-product,  we  apply  the  previously-defined  function  ipn  after 
performing  row  and  column  transforms. 
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Theorem  5.5.1.  Let  N — an,  M = j3m,  with  a and  (3  being  fixed  integers.  Let 

« - mHw{N)mw{M)  (257) 

denote  the  Walsh  transform  of  a two  dimensional  N x M function  f.  Then  g can  be 
expressed  in  terms  of  the  p-product  as: 

g = (W1  ©a  (w2  ©„(...  (w„©„f)))),  (258) 


with 

f = Hw(M)f'  = VK(Vl  ffi„(v2  ©„(...  (vm  ©/')))).  <259) 

where  w,  and  v,  are  the  row  vectors  corresponding  to  F(a)  and  F (/?),  respectively. 
Proof:  Since  the  Walsh  transform  is  separable,  the  proof  is  similar  to  the  one  dimensional 
case. 

Q.E.D. 

The  next  theorem  generalizes  this  result. 


Theorem  5.5.2.  Let  N — p\pi ...  pn,  M = qiq2  . . . qm,  and  f be  an  N x M array.  The 
p-product  factorization  of  the  Walsh  transform  off  is  given  by 

g = A Tm^M  (Wl  (W2  ®P2  C ' ' (Wn  ©p-^))))’  (260) 

with 

f = v>«(vi  ©„  (v2  ©„  (. . . (v,  ©,„f')))),  (261) 


where  w,  and  v,  are  the  row  vectors  corresponding  to  F(pt)  and  F (qf),  respectively. 
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5.6.  /^-product  Matrix  Transforms 


We  have  discussed  the  p-product  formulations  of  the  Fourier  transform,  Walsh 
transform  and  the  generalized  Walsh  transform.  All  of  these  are  orthogonal  transforms. 
Orthogonal  matrices  play  a major  role  in  the  study  of  spectral  analysis  and  in  efficient 
decomposition  schemes  of  unknown  functions  for  computation  on  a digital  computers. 
In  the  conclusion  of  this  chapter,  we  present  /7-product  formulations  of  a class  of 
general  orthogonal  and  nonorthogonal  matrix  transformations  referring  as  p-product 
matrix  transforms. 

n— 1 

For  any  integer  N = II  Pn  we  can  define  a set  of  core  matrices  for  each  factor 

«=o 

Pi.  There  will  be  N/p,  core  matrices  of  dimension  p,  x p,  for  each  factor  p,.  According 
to  Good  [12],  each  core  matrix,  denoted  Ml,  is  defined  as 


M[  = 


m 


m 


i,0,0 


i 

i,p,-l,0 


m 


l 

i,0,  p,  — 1 


m 


*, Pi-1, Pi- 1 J 


(262) 


of  dimension  pi  x pi,  where  / = 0, 1, . . . , ( N/pi ) — 1,  i = 0, 1, . . . , n — 1,  / is  a superscript 
and  not  a power,  and  m\  st  represents  a dummy  variable.  The  core  matrix  will  then  be 
used  to  define  a set  of  n Good  matrices,  so  named  for  their  inventor.  Let  ct(j,  l)  denote  the 


/th-column  of  M|  and  r,  = jr  — 1.  Each  Good  matrix  G;  is  an  N x N matrix  defined  as 


N 


G,  = 


Ci(0,0)  O ...  O Cj(l,0) 

O Ci(0,l)  ...  O O 


O 


o 


Ci(0,r,) 


O 


o 

o 

Ci(i,n) 


...  Ci(pi  — 1,0)  ...  O 

o ...  o 


o 


...  c i(pi  - l,r<). 


(263) 


Note  that  G,  contains  only  Npi  nonzero  entries. 

The  final  general  transformation  matrix  is  defined  as  the  matrix  product  of  the 
respective  Good  matrices  [1];  that  is, 

71  — 1 

h„ = n G« 

t=0 


(264) 
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Specifically,  if  each  G,  is  generated  by  one  core  matrix  M,  (i.e.,  M(  = M™  = M, 
for  all  / and  m),  then  it  has  been  proved  [1]  that 

n— 1 

Hn  = G;  = M0  ® Mi  ® . . . ® Mn_i.  (265) 


i=0 

We  call  H„  a tensor  matrix  product  transform. 

We  now  induce  the  p-product  version  of  this  transform.  We  begin  by  looking  at 
two  examples. 


Example  5.6.1.  If  N = 2 x 2,  then  by  Equation  (264),  the  tensor  matrix  product 
transform  is 


H2  = Mo  ® Mi  = (Mo  <8>  I2XI2  ® Mi) 

= (M0  <g>  I2)P(iV,  2)(Mi  ® I2)P {N,  2). 
We  know  P(Ar,  2)P(Ar,  2)  = I/y.  By  letting 


(266) 


T2  = (M0<g>Mi)P(jV,2)  = H2P(7V,2),  (267) 


it  follows  that 


T2  = (M0®I2)P(JV,2)(Mi®I2).  (268) 

Applying  T2  to  a 4— point  function  f,  we  have  that 

T2f  = (roro(Mo)  ©2(rou;(Mi)  ©2f))',  (269) 

where  roto(M,)  is  the  row  vector  corresponding  to  the  matrix  Mf  . Thus,  we  call  T2f 
the  4 —point  2 -product  matrix  transform  of  f. 
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Example  5.6.2.  Suppose  N = 2x2x2.  By  Equation  (264),  the  tensor  matrix  product 

transform  is 

H3  = Mo  <8>  Mj  <8)  M2 

= (M0  <8>  I4)(l2  <8>  (Mi  ® M2)) 

= (M0  ® I4)P(Ar,  2)((M!  ® M2)  (8)  I2)P(Ar,  4) 

= (Mo  (8)  I4)P(Ar,  2)(Mi  (8)  I4)(P(4, 2)  ® I2)(M2  ® I4)(P(4, 2)  ® I2)P(iV,  4). 

(270) 

It  can  be  proved  that 

(P(4, 2)  ® I2)P(8, 4)(I2  ® P(4, 2))P(8, 2)  = I8.  (271) 


By  letting 


we  have 


T3  = (Mo  ®Mi®  M2)(I2  ® P(4, 2))P(8, 2) 
= H3(I2  <8>  P(4, 2))P(8, 2), 


(272) 


T3  = (M0  (8)  I4)P(8, 2)(M1®I4)(P(4, 2)  ® I2)(M2  ® I4).  (273) 

Applying  this  to  an  8-point  function  f,  we  obtain 

T3f  = (rou>(M0)  02(rou;(M1)  ©2(rou;(M2)  ©2f)))'.  (274) 

We  call  T3f  the  8-point  2-product  matrix  transform  of  f. 

If  we  let  Q be  a permutation  matrix  and  Q = P(8,4)(I2  ® P(4, 2))  in  the  above 
formulation,  then  T3  = H3Q. 

Direct  computation  shows  that  Q satisfies  the  condition 

Q(xi  <8>  x2  ® x3)  = x3  ® x2®xi,  (275) 

where  xi,x2,  and  x3  are  2-dimensional  vectors.  Since  these  tensor  products  span  C8, 
condition  (275)  uniquely  defines  Q.  The  matrix  Q is  called  the  8-point  bit-reversal  matrix. 


77 


In  general,  if  N = 2n,  then  one  can  show  that  the  //-point  bit-reversal  matrix  is 
given  by 

n 

Qw=n(i2-'®p(2i.2))  (276) 

= (I2-2  ®P(4,2))(I2»— 3 ® P(8,2))...  (l2®P(2n-1,2))P(2n,2). 

By  the  induction  on  transform  size,  we  have  the  following  result. 


Theorem  5.6.3.  If  N = 2n,  then  the  N-point  2-product  matrix  transform  is 
Tn  = UnQ(N) 

= (M0®Mi®...®Mn_i)f  n(l2n-.  ®P(2\2)) 

\*=2  > 

= (n  (M,(g)I2n-i)(P(2n-,',2)  ®I2i)j(M„_1®I2»-i). 


(277) 


Proof:  We  know  that  the  theorem  holds  for  n= 1 and  2. 

Assume  that  this  is  also  true  for  k — n — 1.  That  is 
Tn_i  =H„_1Q(2n-1) 

(n— 3 

J][  (M,-  ® I2n-2 ) (P (2n-,_1 , 2)  ® 1 2i)  ) (Mn_2  (8)  I2»-2). 

Using  induction,  we  prove  the  theorem  for  n.  That  is 
Tn  = HnQ(N) 

= (M0®Mi®...®M  (I2n-.  ®P(2‘,2)) 

= (Mo  ® Mi  ® . . . ® M„_i)  ^I2  ® (j[  (l2"-*-i  ® P(2\  2))  j ^ P(2n,  2) 
= (M0  ® Mi  ® . . . ® Mb_!)(I2  (8)  Q(2n"1))P(2n,  2) 

= (Mo  (8)  I2-i)(I2  ®M1®...®M„_i)(I2®  Q(2n“1))P(2n,2) 

- (M0  <8)  I2-i)(l2  ® (Mi  ® . . . ® M„_1)Q(2”-1))P(2n,2) 

= (M0  ® I2r.-i)(I2  <8>  T„_i)P(2n,  2) 

= (M0®I2»-i)P(2b,2)(Tb_1  <8>I2)P(2”,2n-1)P(2n,2) 

= (Mo  ® I2n-1)P(2n,2)(Tn_1  <8>  I2). 


(278) 


(279) 
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Substituting  T„_i  as  given  by  Equation  (278)  results  in  the  desired  conclusion. 


Corollary  1.  Suppose  N — 2"  and  f is  a function  of  length  N,  then  the  2-product  matrix 
transform  of  f is 


T„f  = (rou>(M0)  @2(rou>(Mi)  ©2(. . . ( ©2(row(Mn_i)  02f)))))'.  (280) 


n — 1 

More  generally,  for  N = n pt,  the  permutation  matrix  Q is  the  mixed  radix 

i=0 

analog  of  bit  reversal.  It  is  given  by  the  factorization 


where  n;  = Yl  Pj- 
i= o 

By  induction  on  the  number  of  factors,  the  following  theorem  can  be  obtained. 


n — 1 

Theorem  5.6.4.  Let  N = [7  Pt-  The  N -point p-product  matrix  transform  is  given  by 

i=0 


Q.E.D. 


(281) 


n—i 


(282) 


where  mo  = l,m,  = pop\ . . 


Proof:  If  n=  1,  then  N = pop\  and 

T2  = (M0  ® Mi)P(Ar,p0) 


= (M0  ® IPl)(Ipo  ® Mi)P(7V,p0) 


(283) 


= (M0  ® IPl)P(Af,p0)(Mi  ® IPo), 


which  satisfies  the  formula. 
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Assume  the  formula  is  correct  for  K = pip2 . . . pn-i  and  we  shall  prove  it  for 
N = p0K. 

T „ = (M0  M„_i)  (l„y  0 P(AT/nj,pn_j+i))J  P(N,p0) 

= (M0  0 Mi  0 ...  0 Mn_i)(IPo  0 Q(A'))P(Ar, p0) 

= (MO0I^)(IPO  0 Mi  0 M2  0 . . . 0 M„_i)(IP0  0 Q(K))P{N,p0) 

= (Mo0lir)(IPo  0)  (Mi  0 ...  0 Mn_i)Q(A'))P(Ar, p0) 

= (M0  0lir)P(iV,po)(Hn_iQ(A:)0lJ,o). 

Substituting  Hn_iQ(/i ) by  the  above  formula  for  K provides  the  desired  result. 

Q.E.D. 

n — 1 

Corollary  2.  Suppose  N — n Pi  and  f is  a function  of  length  N,  then  the  p-product 

i= 0 

matrix  transform  of  f is 


Tnf  = (row(Mo)  ©Po  (row(Mi)  ©pi  (. . . ( ©Pn_2  (rou;(M„_i)  ©p^j)))))'.  (285) 


Note  from  the  above  discussion  that  the  only  difference  between  the  tensor  matrix 
product  transform  formulation  and  the  p-product  matrix  transform  formulation  is  the 
permutation  matrix  Q.  But  in  most  signal  processing  transforms,  such  as  the  fast  Fourier 
transform,  the  fast  Walsh  transform  and  the  fast  generalized  Walsh  transform,  the  bit- 
reversal  reordering  processing  is  always  necessary.  Hence,  the  p-product  transform 
expression  provides,  in  general,  for  a more  efficient  computational  framework  than  does 
the  tensor  product  expression. 

Specific  transformations  which  are  readily  implementable  in  the  context  of  the 
above  p-product  factorization  techniques  include  the  Fourier,  the  Walsh,  and  the  general- 
ized Walsh  transforms,  as  discussed  previously,  and  a variety  of  other  similar  transforms. 
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For  example,  if  we  have 


M,  = 


cos  9 
sin  9 


sin  9 
— cos  9 


(286) 


for  all  i,  then  Tn  generated  by  M,  is  an  orthogonal  p-product  matrix  transform.  As  9 
varies  from  0°  to  45°,  it  results  in  different  kinds  of  p-product  matrix  Tn.  In  particular, 
when  9 = 45°,  Tra  is  equivalent  to  the  Walsh  transform  of  order  2. 

In  conclusion,  we  note  that  we  can  use  the  p -product  to  establish  simplified 
algorithms  for  efficient  implementation  of  matrix  operations  over  a class  of  generalized 
tensor  matrix  products.  This  means  that  any  orthogonal  or  nonorthogonal  transforms 
obtainable  from  a series  of  matrix  tensor  products  can  be  expressed  concisely  in  terms 
of  the  p-product  and  may  be  efficiently  implemented  for  computer  applications. 


CHAPTER  6 

THE  GENERALIZED  MATRIX  PRODUCT  AND  THE  WAVELET  TRANSFORM 


The  subject  of  wavelet  analysis  has  recently  drawn  a great  deal  of  attention  from 
mathematical  scientists  in  various  disciplines.  It  is  creating  an  area  of  common  interest 
for  mathematicians,  physicists,  and  electrical  engineers.  We  begin  by  briefly  describing 
the  basic  concepts  of  the  theory  of  wavelet  analysis. 

It  is  well  known  that  the  Fourier  transform  decomposes  a signal  into  individual 
frequency  components  but  does  not  provide  information  as  to  when  the  frequencies 
occurred.  When  the  signal  to  be  analyzed  is  nonstationary,  a relevant  analysis  calls  for 
keeping  the  time  information  in  order  to  exhibit  its  time-varying  spectral  properties.  The 
most  straightforward  solution  is,  therefore,  to  split  the  signal  into  fractions  within  which 
the  stationary  assumptions  apply.  The  Gabor  Transform  (GT)  (or  the  Short  Time  Fourier 
Transform  (STFT))  is  commonly  used  to  perform  this  decomposition.  It  introduces  a 
time-localization  window  function,  g(t  — b),  where  the  parameter  b is  used  to  translate 
the  window  in  order  to  cover  the  whole  time-domain  for  extracting  local  information  of 
the  Fourier  transform  of  the  signal.  The  principal  problem  here  is  that  any  one  choice  of 
g{t)  results  in  windows  that  are  too  wide  to  capture  all  nonstationary  behavior,  and  too 
narrow  to  capture  low-frequency  information.  The  recently  introduced  wavelet  transform 
is  an  alternative  tool  that  deals  with  nonstationary  signals.  The  decomposition  is  carried 
out  by  means  of  a special  analysis  function  if),  called  the  basic  wavelet,  which  is  translated 
in  time  (for  selecting  the  part  of  the  signal  to  be  analyzed),  then  dilated  or  contracted 
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using  a scale  parameter  (in  order  to  focus  on  a given  range  of  oscillations).  In  addition 
to  inheriting  the  basic  properties  of  g,  ip  satisfies  an  additional  condition 

OO 

J ip(x)dx  — 0.  (287) 

— OO 

This  latter  condition  provides  an  extra  degree  of  freedom  for  introducing  a dilation  (or 
scale)  parameter  in  order  to  make  the  time-frequency  window  flexible.  Analogous  to 
Fourier  analysis,  corresponding  to  the  integral  wavelet  transform,  there  is  a wavelet 
series.  While  the  integral  wavelet  transform  is  defined  to  be  the  convolution  with  respect 
to  the  dilation  of  the  reflection  of  some  function  ip,  the  wavelet  series  is  expressed  in 
terms  of  a single  function  ip,  called  an  R-wavelet  (or  simply,  a wavelet)  by  means  of 
two  very  simple  operations:  binary  dilations  and  integral  translations.  However,  unlike 
Fourier  analysis,  the  integral  wavelet  transform  with  a basic  wavelet  ip  and  the  wavelet 
series  in  terms  of  a wavelet  ip  are  intimately  related.  In  fact,  if  ip  is  chosen  to  be  the  dual 
of  ip,  then  the  coefficients  of  the  wavelet  series  of  any  square-integrable  function  / are 
precisely  the  values  of  the  integral  wavelet  transform,  evaluated  at  the  dyadic  positions  in 
the  corresponding  binary  dilated  scale  levels.  This  definition  provides  the  integral  wavelet 
transform  of  a function  / to  simultaneously  localize  / and  its  Fourier  transform  / with 
zoom-in  and  zoom-out  capability.  There  also  exists  real-time  algorithms  for  obtaining  the 
coefficient  sequences  of  the  wavelet  series,  and  for  recovering /from  these  sequences.  It 
is  for  this  reason  that  this  new  transform  has  drawn  a great  deal  of  attention,  in  a very 
short  time  span,  from  mathematicians  and  scientists  in  various  disciplines.  Its  numerous 
applications  can  be  found  in  [3,  6,  16,  17]. 

Recall  from  the  previous  chapters  that  the  p-product  provides  novel  algorithms 
for  computing  and  expressing  Fourier  and  Walsh  transforms.  In  this  chapter,  we  will 
introduce  the  p-product  into  the  formulation  of  the  wavelet  transform.  A well  known  fact 


83 


of  linear  algebra  is  that  linear  transforms  can  be  represented  in  terms  of  matrix-vector 
products.  However,  since  the  wavelet  transform  is  a function  of  two  variables  (both 
time  and  frequency),  it  is  difficult  to  express  this  transform  in  matrix  product  form.  Even 
though  Heller  et  al.  [6]  have  defined  wavelet  matrices,  they  are  only  used  to  prove  certain 
mathematical  properties  but  not  as  a transform  representation.  In  contrast,  the  generalized 
matrix  product  provides  for  expressing  the  wavelet  transform  in  matrix  form  succinctly. 
In  this  chapter,  we  show  how  to  represent  the  wavelet  transform  and  its  inverse  in  terms 
of  the  p-product  by  using  the  wavelet  matrices  defined  in  [6].  In  addition,  we  describe  a 
simple  and  fast  wavelet  transform  algorithm  using  the  p-product  and  parallelism. 

After  providing  a brief  introduction  of  the  wavelet  transform,  a new  algorithm  for 
computing  the  wavelet  transform  in  terms  of  the  p-product  is  established.  The  algorithm 
is  based  on  the  decomposition  of  a long  summation  into  several  short  ones  and  then 
uses  the  p-product  to  carry  out  the  computation.  This  provides  a fast  algorithm  for 
wavelet  transform  computation.  When  executing  this  algorithm  on  parallel  machines, 
the  computing  time  is  at  least  g times  less  than  the  standard  method,  where  g denotes 
the  genus  of  the  wavelet  matrix.  Predictively,  the  efficiency  improves  following  the 
architecture  of  parallel  machines  employed. 

6.1.  The  Wavelet  Transform 

In  this  section,  we  give  a brief  introduction  of  the  wavelet  transform  and  define 
some  basic  notations  that  are  going  to  be  used  throughout  this  chapter.  The  details  can 
be  found  in  [5]. 
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6.1.1.  From  Fourier  Analysis  to  Wavelet  Analysis 

Let  L2(0,27r)  denote  the  set  of  all  measurable  27r -periodic  functions  of  / defined 
on  the  interval  (0,  2x)  with 


2x 

J \f(x)\2dx  < oo. 

0 

Any  function  /in  L2(0,27r)  has  a Fourier  series  representation 


(288) 


/(*)=  J]  c„einx , 

n=— oo 

where  the  coefficients  cn  are  defined  by 

2 TT 

c’  = h/nx) 


(289) 


e~tnxdx. 


(290) 


The  convergence  of  the  series  in  Equation  (289)  is  in  L2(0,27r),  namely: 


2J  N 

lim  \f(x)~  cnetnx\2dx  = 0. 

M,N—>oc  J L — ' 


(291) 


n——M 


It  follows  that / can  be  decomposed  into  a sum  of  infinitely  many  orthogonal  components 
gn(x)  :=  cne,ni,  where  orthogonality  is  defined  by 

2x 

{gm,  gn)  = ^ J 9m(x)gn(x)  dx  = 0,  for  all  m±n.  (292) 


It  is  obvious  that 


wn[x)  - etnx,  n = 1,0,1,...  (293) 

is  an  orthonormal  (o.n.)  basis  for  L2(0,  27t),  which  is  generated  by  the  dilates  of  a single 
function 


w(x)  = e,x; 


(294) 
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that  is,  wn{x)  = w(nx)  for  all  integers  n.  This  means  that  every  27r -periodic  square- 
integral  function  is  generated  by  a superposition  of  integral  dilations  of  the  basic  function 
w(x)  — elx.  Hence,  we  call  the  basic  function 

w(x ) = e,x  = cos  x + i sin  x (295) 

a sinusoidal  wave.  Note  that  this  is  the  only  function  required  to  generate  all  2t- 
periodic  square- summable  functions.  For  any  integer  n with  large  absolute  value,  the 
wave  wn(:r)  = w{nx ) has  high  frequency,  and  for  n with  small  absolute  value,  the  wave 
wn  has  low-frequency.  Therefore,  any  function  in  L2(0,27r)  is  composed  of  waves  with 
various  frequencies. 

We  next  consider  the  space  L2(R)  of  measurable  functions  /,  defined  on  the  real 
line  R,  that  satisfy 

OO 

J | f(x)\2dx  < oo.  (296) 

— OO 

Clearly,  the  two  function  spaces  L2(0,27 r)  and  L2(R)  are  quite  different.  In  particular, 
since  (the  local  average  value  of ) every  function  in  L2(R)  must  decay  to  zero  at  +oo,  the 
sinusoidal  (wave)  functions  wn  do  not  belong  to  L2(R).  In  fact,  if  we  look  for  waves  that 
generate  L2(R),  these  waves  should  decay  to  zero  at  Too,  and  for  all  practical  purposes, 
the  decay  should  be  very  fast.  That  is,  we  look  for  small  waves,  or  wavelets,  to  generate 
L2(R).  As  in  the  situation  of  L2(0,27t),  where  single  function  w(x)  = etx  generates 
the  entire  space,  in  L2(R),  we  use  a set  of  basis  functions  which  are  discrete  scales  and 
translates  of  a single  basis  wavelet  ip  : 


fpj,k(x)  = a02/2V> 


, j,  k G Z,a0  > Mo  7^  0. 


(297) 
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In  recently  published  papers,  authors  prefer  to  use  a dyadic  discrete  wavelet  scheme 
(i.e.,  ao  — 2).  We  follow  this  trend  in  this  section  by  employing  integral  powers  of  2 for 
frequency  partitioning.  As  a consequence,  we  consider  small  waves 


(i.e.,  dilation  by  23)  and  a dyadic  translation  (of  k/23).  Thus,  we  are  interested  in  wavelet 
functions  ip  whose  binary  dilations  and  dyadic  translations  are  enough  to  represent  all 
the  functions  in  L2(R). 

We  begin  by  considering  an  orthogonal  basis  generated  by  ip.  In  the  following  the 
symbol  6tJ,  where  i,j  £ Z,  denotes  the  6-function;  i.e., 


Definition  6. 1.1.1.  A function  ip  £ L2(R)  is  called  an  orthogonal  wavelet  (or  o.n. 
wavelet),  if  the  family  {fpj,k}>  defined  by 


ip(23x  — k) , j,  k 6 Z. 


(298) 


Note  that  ip(23x  — k ) is  obtained  from  a single  wavelet  function  ip(x)  by  a binary  dilation 


(299) 


(300) 


is  an  orthonormal  basis  of  L2(R);  that  is, 


(,^j,kt  — 67  / ■ bk,mi  ji  k,l,T7l  £ Z, 


(301) 


and  every  f £ L2(R)  can  be  written  as 


OO 


(302) 
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where  the  convergence  of  the  series  in  Equation  (302)  is  in  L2(R),  namely: 


lim 

Ml  ,N\M2,N2—KX> 


n2 

Ni 

f-  E 

y.  ch^xi)j,k 

j=—M2  k=—M\ 

= 0. 


(303) 


The  simplest  example  of  an  orthogonal  wavelet  is  the  Haar  function  tpjj  defined  by 


( 1 for  0 < x < 1/2, 

*I>h(x)  ■=  \ -1  for  1/2  < x < 1,  (304) 

1 0 otherwise. 


The  series  representation  of  / in  Equation  (302)  is  called  a wavelet  series. 
Analogous  to  the  notion  of  Fourier  coefficients  in  Equation  (290),  the  wavelet  coefficients 
Cj^  are  given  by 

Cj,k  = (f^j,k)-  (305) 

That  is,  if  we  define  an  integral  transform  on  L2(R)  by 

OO  

(W^f)(b,a):=\a\~i  J f € L2(R),  (306) 

— OO 

then  the  wavelet  coefficients  in  Equations  (302)  and  (305)  become 

•w-Pw)  (*■*)•  <307> 

The  linear  transformation  is  called  the  integral  wavelet  transform  relative  to  the  basic 
wavelet  ip.  Hence,  the  (j,  k)th  wavelet  coefficient  of  / is  given  by  the  integral  wavelet 
transformation  of /evaluated  at  the  dyadic  position  b = k/2]  with  binary  dilation  a = 2-J  , 
where  the  same  o.n.  wavelet  ip  is  used  to  generate  the  wavelet  series  (302)  and  to  define 
the  integral  wavelet  transform  (306). 
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Notice,  first,  that  this  integral  transform  greatly  enhances  the  value  of  the  (integral) 
Fourier  transform  T,  defined  by 

OO 

(Ff)(y)  :=  J e~ixy f(x)dx,  f € L2(R).  (308) 

— OO 

As  is  well  known,  the  Fourier  transform  is  the  other  important  component  of  Fourier 
analysis.  Hence,  it  is  interesting  to  note  that  while  the  two  components  of  Fourier 
analysis,  namely:  the  Fourier  series  and  the  Fourier  transform,  are  basically  unrelated; 
the  two  corresponding  components  of  wavelet  analysis,  namely:  the  wavelet  series  (302) 
and  the  integral  wavelet  transform  (306),  have  an  intimate  relationship  as  described  by 
Equation  (307).  Second,  the  integral  wavelet  transform  (W^f)(b,a)  gives  the  location, 
the  rate  (in  terms  of  a),  and  the  amount  (measured  by  the  value  (VF^/)(6,  a))  of  change 
of  /,  with  the  zoom-in  and  zoom-out  capability.  This  means  that  the  wavelet  transform 
gives  a time-scale  representation  which  is  sharp  in  time  at  small  scales  and  sharp  in 
scale  at  large  scales.  This  information  is  extremely  valuable  in  many  applications  such 
as  time-frequency  analysis.  For  instance,  in  data  compression,  the  values  of  (W^f)(b,  a) 
below  a certain  tolerance  level  are  removed,  while  in  lowpass  filtering,  (W^/)(&,  a)  is 
replaced  by  zero  for  small  values  of  a. 

Analogous  to  the  Fourier  transform,  the  function  /can  be  reconstructed  from  the 
values  of  (W^f)(b,a).  Any  formula  that  expresses  every  / e L2(R)  in  terms  of 
(W^f)(b,a)  will  be  called  an  inverse  formula.  There  are  several  forms  of  this  inverse, 
which  can  be  found  in  [5].  Here,  we  only  give  the  most  general  inverse: 

OO  OO 

s(x)  = h / / <309) 

— oo  0 

The  (kernel)  function  to  be  used  in  the  formula  of  (W^f)(b,a)  will  be  called  a dual 
of  the  basic  wavelet  0.  Hence,  in  particular,  0 can  be  used  as  a basic  wavelet,  only  if 


an  inversion  formula  exits. 
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6.1.2.  The  Integral  Wavelet  Transform  and  Time-Frequency  Analysis 

It  is  well  known  that  the  Fourier  transform  T is  not  only  a very  powerful 
mathematical  tool,  but  also  has  physical  significance.  For  instance,  if  a function 
/ e L2(R)  is  considered  as  an  analog  signal  with  finite  energy,  defined  by  its  norm 
||/||2,  then  the  Fourier  transform 

f(u)  :=  (IF /)(u>)  (310) 


of  / represents  the  spectrum  of  this  signal.  In  signal  analysis,  analog  signals  are 
defined  in  the  time-domain,  and  the  spectral  information  of  these  signals  is  given  in 
the  frequency-domain.  If  we  allow  negative  frequencies,  then  both  domains  are  the  real 
line  R. 

However,  the  formula 

OO 

f(u)=  J e~itu} f(t)dt  (311) 

— OO 

of  the  Fourier  transform  alone  is  quite  inadequate  for  most  applications.  In  the  first 
place,  to  extract  the  spectral  information  f(u)  from  the  analog  signal  f(t ) by  using  this 
formula,  it  takes  an  infinite  amount  of  time,  using  both  past  and  future  information  of  the 
signal  just  to  evaluate  the  spectrum  at  a single  frequency  u.  Besides,  the  formula  does 
not  even  reflect  frequencies  that  evolve  with  time.  What  is  really  needed  is  the  ability  to 
determine  the  time  intervals  that  yield  the  spectral  information  on  any  desirable  range  of 
frequencies  (or  frequency  band).  In  1946,  Gabor  developed  a frequency  decomposition, 
which  is  also  referred  to  the  Short-Time-F ourier  transform  or  Gabor  transform.  He  used 
a Gaussian  function 


9a(t)  := 


5 


(312) 
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where  a> 0 is  fixed,  as  the  window  function  to  give  a time-localization  of  the  Fourier 
transform.  For  any  fixed  value  of  a>0,  the  Gabor  transform  of  an  / £ L2(R)  is  defined 
by 

OO 

($"/)(«)=  / (e““'/(i))9«(‘  - b)dt.  (313) 

— OO 

That  is,  (Q% /)  (u)  localizes  the  Fourier  transform  of  / around  t=b.  The  width  of  the 
window  is  determined  by  the  (fixed)  positive  constant  a to  be  discussed  later.  Observe 

OO  OO 

that  taking  u=0  and  a = (4a)-1,  we  have  / ga[t-b)db=  J ga(x)dx  = 1.  Thus, 

— OO  — OO 

OO 

J (etf)(u)db  = f( CJ),  »eR.  (314) 

— OO 

That  is,  the  set  {G£f  : b e R}  of  Gabor  transforms  of  /decomposes  the  Fourier 
transform  f of  f exactly,  to  give  its  local  spectral  information. 

However,  since  the  frequency  of  a signal  is  directly  proportional  to  the  length  of 
its  cycle,  it  follows  that  for  high-frequency  spectral  information,  the  time-interval  should 
be  relatively  small  to  give  better  accuracy,  and  for  low-frequency  spectral  information, 
the  time-interval  should  be  relatively  wide  to  give  complete  information.  In  other  words, 
it  is  important  to  have  a flexible  time-frequency  window  that  automatically  narrows  at 
high  center-frequency  and  widens  at  low  center-frequency.  As  pointed  out  above,  the 
Gabor  transform  yields  a time-frequency  decomposition  with  a fixed  bandwidth.  Thus,  it 
cannot  provide  a satisfactory  performance  in  this  case.  Fortunately,  the  integral  wavelet 
transform  relative  to  some  basic  wavelet  ip,  as  introduced  above,  has  this  so-called 
zoom-in  and  zoom-out  capability. 

To  be  more  specific,  both  ip  and  its  Fourier  transform  ip  must  have  sufficiently 
fast  decay  so  that  they  can  be  used  as  window  functions.  For  an  L2(R)  function  w to 
qualify  as  a window  function,  it  must  be  possible  to  identify  its  center  and  width,  which 


are  defined  as  follows. 
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Definition  6. 1.2.1.  A nontrivial  function  we  L2(R)  is  called  a window  function  ifxw(x) 
is  also  in  L2(R).  The  center  t*  and  radius  Aw  of  a window  function  w are  defined  to  be 


OO 

**  :=  7^2  / xKx)r 

No  J 


dx 


(315) 


and 


1/2 


(316) 


Aw  Pt  | / (x~rf\wW\2dx 

\— OO 

respectively,  while  the  width  of  the  window  function  w is  defined  by  2AW. 

Note  that  Gabor  transform  is  a special  case  of  wavelet  transforms.  Applying  the 


above  definition,  the  width  of  the  Gabor  transform  window  function  ga  is 

1/2 


2A‘“ = M 


J x2gl(x)dx 


= 2 \fa. 


(317) 


6.2.  Wavelet  Matrices 

Following  the  introduction  of  the  wavelet  transform  in  last  section,  we  now  come 
to  define  the  finite  wavelet  transform  by  using  the  wavelet  matrices  defined  in  [6]. 

In  1988,  Ingrid  Daubechies  defined  the  notion  of  the  multiplier  2 compactly  sup- 
ported discrete  wavelet  transform  and  obtained  conditions  for  smoothness  and  polynomial 
representation  by  multiplier  2 wavelet  series  [10].  In  particular,  she  defined  a scaling 
function  <f>(x)  as  a compactly  supported  solution  of 

2g-l 

<KX)  = XI  akH2x  - k)i  (318) 

*=0 

where  a0,  a1} . . . , a2g- 1 are  the  scaling  coefficients.  Associated  with  this  scaling  function 
of  the  wavelet  system  there  is  another  set  of  coefficients,  6*  = (— l)ka2g-i-k,  which 
defines  the  wavelet  function  as 

29-1 

^(x)  = X - *)• 

h=0 


(319) 
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Using  these  definitions,  Heller  et  al.  [6]  introduced  wavelet  matrices  as 
generalizations  of  the  2x2 g matrix  of  the  form 


where  the  a's  and  b's  are  defined  as  above.  It  is  not  difficult  to  ascertain  that  this  satisfies 
the  wavelet  scaling  property  and  that 


where  m is  the  rank  of  the  matrix,  and  g denotes  the  genus  of  the  wavelet  matrix;  i.e., 
the  number  of  m x m blocks  in  the  matrix.  The  vector  a0  = a®, . . . , a^_1^  is 

called  the  scaling  vector  and  for  0 < s < m,  the  5th  row  of  a,  denoted  by  as,  is  called 
a wavelet  vector. 

Note  that  wavelet  matrices  of  rank  m will  correspond  to  wavelet  systems  with 
multiplier  m,  replacing  the  multiplier  2 used  by  both  Daubechies  and  Mallat  in  [10,  18]. 


(320) 


(321) 


it 


k 


In  order  to  generalize  this  concept  one  may  define 


(322) 


where  a-1  = a,  and  a]  = bt.  The  general  m x mg  matrix  is  then  of  form 

/ «8  «?  •••  <,-i\ 

«0  G1  • • • amg- 1 


(323) 


with  the  wavelet  scaling  conditions 


k 


(324) 


k 
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Analogous  to  the  case  of  m= 2,  under  this  definition  of  the  wavelet  matrix  a,  a scaling 
function  0°  and  m- 1 fundamental  wavelets  are  defined  by  the  system  of  equations 

ipr(x)  = akip°(mx  — k ),  (325) 

it 

where  0 < r < m.  The  scaling  function  is  to  be  thought  of  as  a low-pass  function 
while  the  fundamental  wavelets  are  high-pass  functions.  If  we  define  a set  of  auxiliary 
functions  by  the  formula 

rpjj,  (mJ x — k ),  (326) 

where  j,  k e Z,  then  Resnikoff  [20]  has  shown  that  the  set 

{*pjk(x)  : 0 < r < m,j,k  E Z}  (327) 

is  an  orthonormal  basis  for  L2(R).  The  support  of  ipTjk(x)  has  length  equal  to  the  length 
of  the  support  of  the  scaling  function  tp°(x)  divided  by  rn] , where  the  quantity  j is  called 
the  scale  of  the  wavelet  function.  This  gives  a corresponding  relationship  between  a 
wavelet  matrix  a and  an  orthonormal  wavelet  basis.  More  details  of  this  discussion  can 
be  found  in  Section  4 of  this  chapter. 

A wavelet  matrix  for  which  g=l;  i.e.,  a square  wavelet  matrix,  is  said  to  be  a Haar 
matrix.  A number  of  classical  examples  of  specific  matrices  which  have  different  origins 
in  mathematics  and  signal  processing  can  all  be  seen  to  be  Haar  wavelet  matrices  of 
specific  types.  These  include  the  finite  Fourier  transform  matrices,  the  discrete  cosine 
transform  matrix,  Hadamard  and  Walsh  matrices,  Rademacher  matrices,  and  Chebyshev 
matrices. 

Let  f(x)  = ^2  fn(x),  where  /„  is  a sequence  of  functions  x — » fn(x)  defined  on 

n 

some  infinite  set  (e.g.,  Z or  R).  The  function  f(x)  will  have  a meaning  that  is  prescribed 
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by  the  type  of  convergence  that  is  assumed.  Let  us  suppose  that  for  each  x,  only  finitely 
many  of  the  numbers  fn(x)  are  nonzero  and  assume  ask  = 0 unless  0 < k < mg.  Heller 
et  al.  [6]  have  proved  the  following  theorem. 

Theorem  6.2. 1.1.  Let  f : Z — > C be  an  arbitrary  function  defined  on  the  integers  and  let 
a be  a compact  wavelet  matrix  of  rank  m and  genus  g defined  by  Equation  (323).  Then 
f has  a unique  wavelet  matrix  expansion 

m— 1 

f(n ) = clan-ml  + ckan-mki  (328) 

l£l  kel  5—1 

where 

c/  = — V /(n)a„_m/,  (329) 

m ' 

n 

and 

4 = 4 Y,  (330) 

n 

The  wavelet  matrix  expansion  is  locally  finite;  i.e.,for  given  n only  finitely  many  terms  of 
the  series  are  different  from  zero. 

Note  that  using  the  language  of  signal  processing,  the  first  term  in  Equation  (328) 
is  the  low-pass  part  of  the  expansion,  and  the  second  term  is  the  high-pass  part  of  the 
expansion. 

6.3.  The  Generalized  Matrix  Product  and  Wavelet  Transforms 

We  now  know  that  corresponding  to  a given  orthonormal  wavelet  basis,  there  exists 
a wavelet  matrix.  The  algebraic  and  geometric  properties  of  the  set  of  wavelet  matrices 
of  finite  rank  have  been  given  in  [6].  In  this  section,  we  are  going  to  employ  this  wavelet 
matrix  to  induce  the  generalized  matrix  product  form  for  wavelet  transforms. 
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According  to  Theorem  6.2. 1.1,  for  any  function  / : Z — > C,  there  is  a unique 
wavelet  matrix  expansion 

oo  m — 1 oo 


where 


f(n)  = S Cla"-ml  ckan-mk > 

/=— oo  s=l  k=—oo 

cl  = ^ ^ f(n)an— mli 

m 

n 


and  a is  a compact  wavelet  matrix  of  rank  m and  genus  g of  the  form 


( «0 

ai 

• amg—  1 

b0 

b\  . 

A 1 

°mg- 1 

w_1 

b™~1  . 

im- 1 
°mg- 1 

(331) 


(332) 


(333) 


Recall  that  the  wavelet  matrix  expansion  is  locally  finite,  we  denote  f = 
[/(0),/(l),...,/(iV  — 1)]  and  assume  ask  = 0 unless  0 < k < mg  while  denoting 
a°k  = ak. 

Additionally,  if  we  denote  c”  = q for  / € Z,  then  two  equations  in  Equation  (332) 
can  be  combined  to  the  following  expression: 


4 = ~ /(nK-mJb>  o < 5 < m.  (334) 

n 

Using  the  above  assumption  on  a,  we  have 

mg+mk—  1 

4 = - E <335> 

m 

n=mk 

Since  theoretically  k ranges  from  — oo  to  +oo,  in  what  follows,  we  analyze  Equation 
(335)  in  terms  of  two  cases;  i.e.,  k > 0 and  k < 0,  and  denote  the  results  as  and 
ci,  respectively. 


96 


For  k > 0,  we  have 


mg-\-mk— 1 


4 


= - E /(")<-. 

n—mk 

= — 52  f(U  + 


(336) 


n=0 


Decomposing  the  summations  in  Equation  (336)  into  g smaller  ones  of  length  m,  we 
obtain 

mg—l 

4 = ~ f(n  + km)K 


m 

n=0 

m — 1 ..  m— 1 


^ 711  — 1 i 711  — 1 

= — y"  f{p  + km)as  + — Y]  f(m  + p + km)asm+  + 
p— o p= o 

^ m— l 

+ •••  + — S f(m(9  ~ 1)  + P + km)asm{g_lHp. 


(337) 


p= o 

Since  the  length  of  f is  N and  f(n ) = 0 for  n > iV,  the  length  of  the  vector  cg_  is 
That  is,  c+  is  of  the  form  eg.  = eg,  cf, . . . c^m_1  . For  0 < j < g,  let 

hj  = [/(jm),  • • • , f(jm  + N -1)]  ©r 


umj 


- am(j+l)-l  . 


(338) 


By  definition  of  the  /7-product,  the  size  of  h;  is  1 x It  is  not  difficult  to  prove  that  the 
first  summation  of  csk  in  Equation  (337)  corresponds  to  ho(k),  the  second  one  to  hi(k), 
and  the  jth  one  to  hj_\(k). 

Now,  we  denote 


Wn  = 


«0 

1 

S>l 

• J." 

i 

• 

, . . . , Wj  — 

-aTO-l  - 

—s 


■»(»— i) 


amg- 1 J 


(339) 


and  let 


h = [/(mfc),  f(mk  + 1),  • • • , f(N  - 1),  0, . . . , 0]  be  of  size  lxiV, 


(340) 
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where  the  support  of  f*  has  length  equal  to  the  length  of  the  support  of  f minus  mk.  This 
means  that  f*.  can  be  obtained  by  shifting  f to  the  left  mk  elements.  Thus, 

hy  - ij  ©mw*.  (341) 


Combining  Equations  (337)  and  (341)  yields: 
g- 1 


J * * l 

c+  = — h;  = ~ [f  + • • • + fff-l  ©mWj-l]  , 0 < s < m. 

m ' m 


(342) 


j=0 


In  the  case  of  k < 0,  we  substitute  k in  Equation  (335)  by  -k,  while  still  assuming  k> 0. 
Then,  Equation  (335)  becomes 

mg—mk— 1 mg—  1 


c~h=m  ^ f(n)K+mk  = —^2f(n~  mk)asn 
n=—mk  n=0 

Since  /(n)  = 0 unless  0 < n < N,  we  have 

mg— 1 


(343) 


-*  = — f(n~  mk 


(344) 


n=mk 


Thus  the  size  of  cs  = 


"-0+1’  * * • ’ — 2’  C-1 


is  1 x (g  — 1),  and  only  the  first  m(g  — 1) 


terms  of  f are  used.  Therefore,  we  denote 


e=[/(0),/(l),...,/(m(p-l)-l)],  (345) 

and 

e_,  = [0, 0, . . . , 0,  /( 0), . . . , f(m(g  - i - 1)  - 1)],  (346) 

which  has  the  same  size  as  e for  1 < i < g — 2.  Analogous  to  f*,  the  support  of  e_, 
has  length  equal  to  the  length  of  the  support  of  e minus  mi.  This  means  that  e_,  can  be 
obtained  by  shifting  e to  the  right  mi  elements.  Following  the  same  procedure  as  above. 
Equation  (344)  can  be  decomposed  as 
cs-  = [cL,+1,...cL2,ci1] 

= — [e  ©mWs— 1 + e-l  ©mW3- 2 + • • • + e~g+ 2 ©mWl]  ) 


(347) 
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for  0 < s < m.  Combining  Equations  (342)  and  (347),  the  wavelet  transform  of  f under 
the  wavelet  matrix  a is  given  by: 


c=[c°,c1>...lc— »]', 


(348) 


where 


c*  = [ci;<4]  = 


rs  r3  r3  • r3  r3 

c-g- f-1  ? • • • * c— 2)  C—1  j c0  ? • * * * 


— ~ [e  ©mwj-l  + • • • + e-g+2  ®mW?;f  ©mw0  + • • • + fff-l  ®mws-l]  > 


(349) 


for  0 < s < m. 


Before  deriving  the  discrete  inverse  wavelet  transform  of  f,  we  prove  the  following 
lemma. 


Lemma  6.3.1.  Let  w,  be  defined  as  Equation  (339).  Then 


e eh®  Mr)  = p-  i r/i 

0<s<m0<i,j<g  ^ 7 J 


(350) 


Proof: 

E E H®K*)') 

0<s<m  0 <i,j<g 


( 

'ami 

\ 

E E 

\ 

® [Gm;  i amj+l » • • • » l] 

0<s<ra  i,j 

V 

_am(i+l)-l  . 

/ 

(351) 


= EE 


UjmiUjmj 


amiamj+ 1 


-am(«+l)-lamj  am(i+l)-lamj+l 


amiam(j+l)—l 


am(i+l)— laro(j+l)—  1 - 


By  the  wavelet  scaling  conditions  given  in  Equation  (324),  we  have 

0<s<m  /gZ 


(352) 


99 


The  result  now  follows. 

Q.E.D. 

We  employ  the  preceding  lemma  in  order  to  derive  the  p-product  wavelet  transform. 

Theorem  6.3.2.  Let  f be  an  N-length  function  defined  on  the  integers  and  let  a be 
a compact  wavelet  matrix  of  rank  m and  genus  g defined  by  Equation  (323).  Let 
e = [/( 0),  /( 1), . . . f(m(g  - 1)  - 1)].  The  wavelet  transform  off  is  given  by 


c=  [c°lc1,...,cro-1],> 


(353) 


where 


c*  = = 

m 


rs  rs  rs  ■ rs  rs  rs 

c—g- (-1 1 • • • > c— 2i  c— 1 > c0  i ci  ) • • • i 


— ~ [e  ®mwg-l  + • • • + e— 9+2  ®mwl'>  f ®mw0  + • • • + ^9-1  ©mWff-l]  i 
for  0 < s < m,  and  its  inverse  transform  is  given  by 

f=  £ d‘. 


(354) 


0<s<m 

where  ds  = Cq  <g)  (w80)'  + c[  <g>  (w *)'  + . . . + csg_l  <g>  (w*_1)/,  (355) 

3 S S 

C—ii  c— i+1  * • • • » c-i+N/m-l 
and  wf  f,-,  and  e_;  are  as  defined  by  Equations  (339),  (340),  and  (346),  respectively. 
Proof:  On  the  above  discussion,  we  have  shown  that  Equation  (354)  is  valid.  Now  we 
prove  that  Equation  (355)  is  true.  Write 


f,= 


v°  V1  V7 

Vl  5 Vl  '>•••')*  I 


(356) 


with  VtJ  = [ fi(mj ), fi{m(j  + 1)  - 1)],  and 


e_,  = 


v°  V 1 vy~ 

v — n v — n * * * 5 y — i 


<7-2 


(357) 


100 


Then  fo  = f,  and  eo  = e.  Following  the  notation  of  f,,  we  have 

V/  = [ fi(mj ), . . . , + 1)  - 1)] 

= + *)),  • • • , f(m(j  + * + 1)  - 1)] 


(358) 


= vj* 


We  know  that 
c = [c_;c+J 

1 

m 


~ — [e  ©mwg-l  + • • • + e_9+ 2 ©mw{ ; f ®mWfl  + • • • + fj-l  ©mw3-l]  • 

(359) 


Using  Lemma  2.1.10  and  substituting  f,  and  e_,  by  V/ , the  above  equation  can  be 
written  as 


cs  = 


1 


m 


■fl-2 

9-2 

9-1 

9-1  AT  , 

^-1-,, 

E • 

O 

.11 

i=0 

fc=0 

Jt=0 

c-ff+l)  c-</+2)  • • 

53  3 

> c0*  C1  > • • • > cK—i 

By  the  definition  of  c?,  we  have 

1 


S 3 3 

C0ic cN/m—l 


m 

rs  — — Ls  rs  rs 

C1  “ m [C-1)C0>  • • • >CJV/m-2 


(360) 


Cf  = — 
m 


C~i+ 1 > • • • > CN/m-i+l 


(361) 


1 r 


S-i  = 


m 


c-g+ 1 ’ c-g+2i  • ■ • > cN/m-g 


Substituting  these  into  da,  we  obtain 

ds  = cj  ® (wg)'  + c[  ® (wj)'  + • • ■ + c’_!  ® (wj_:)' 
r?-1  9- 1 


1_ 

m 


+ ...+ 
m 


E (w*  ® m')  * • • • * E _1  (w*  ® (wS)') 

U=0  k—0 

e H-i-.  ® K-i)') ^ H ® (w;-i)0 

L»=o 


m 


k=0 


(362) 


9-1 

9-1 

<7  — 1 

£ if-j 

(wj  ® (W;)')  , E (w?  ® (wj)')  , 

.*,i=0 

*,9=0 

*,9=0 
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and 


(363) 


By  Lemma  6.3.1,  we  have  that: 


(364) 


S 


Q.E.D. 


There  are  several  important  observations  we  need  to  make. 

First,  note  that  the  expressions  of  the  wavelet  transform  of  f and  its  inverse  are 
simple  and,  with  the  exception  of  the  p-product,  the  only  operations  required  are  shifts, 
which  can  be  easily  accomplished  on  any  digital  machine.  In  Equation  (341),  there  are 
g p-products  and  each  p-product  requires  N operations,  while  in  Equation  (347),  there 
are  g — 1 p-products  and  each  p-product  requires  m(g  — 1)  operations.  Hence,  the  total 
number  of  operations  required  for  computing  cs  is  Ng  + m(g  - 1)  . However,  since  each 
p- product  is  independent  of  the  others,  after  shifting  f to  form  fi, . . . , fg_i,  and  shifting  e 
to  form  e_i,  e_2, . . . , e~g+2,  all  p-product  computations  can  be  executed  simultaneously 
in  parallel,  as  shown  in  Figure  5.  Therefore,  the  computing  time  of  the  p-product  wavelet 
transform  depends  only  on  the  length  of  the  function  f,  which  is  at  least  g times  faster 
than  the  existing  methods  that  exhibit  complexity  0(mgN)  [3,  6,  5,  10].  Furthermore, 
since  the  support  of  f,  (or  e_,)  equals  the  length  of  the  support  of  f (or  e)  minus  mi,  the 
exact  computing  time  will  be  less  than  predicted  by  theory. 

Second,  the  size  of  cf_  is  1 x (g  — 1),  and  it  affects  only  the  left  boundary  of  the 
transform  cs.  If  g and  m are  very  small,  which  happens  in  most  applications,  or  the 
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boundary  effect  is  less  important,  which  is  the  case  on  big  digital  images,  then  we  can 
ignore  cf_  and  set  it  to  0.  It  follows  that  fewer  processors  and  spaces  will  be  needed  and 
the  formula  will  become  simpler  [30]. 

Third,  the  data  arrangement  of  Equation  (365)  provides  an  advantage  of  implement- 
ing the  wavelet  transform  on  hypercube  parallel  machines.  The  arrangement  is  similar  to 
the  FFT.  Since  the  size  of  each  w*  is  m x 1,  by  the  parallel  property  (Theorem  1.2.10)  of 
the  p-product,  if  we  execute  each  m-product  f^  @mw \ in  parallel,  as  few  as  m operations 
are  required.  Thus,  the  complexity  of  this  new  transform  is  independent  of  the  length  of 
f.  This  is  extremely  valuable  for  large  data  signal  process.  The  larger  the  size  of  f,  the 
more  time  can  be  saved.  However,  there  is  trade-off.  The  number  of  processors  required 
depends  on  the  size  of  f,  that  is,  upon  ^ . Thus,  the  implementation  and  efficiency  of  this 
new  wavelet  transform  formulation  rely  heavily  upon  the  target  architecture.  In  general, 
we  can  save  time  by  at  least  a factor  g.  With  sufficient  parallelism,  the  efficiency 
increases  to  m vs  0(mgN)  [6]. 


Figure  5.  Parallel  computation  of  c+  from  Equation  (342). 

Finally,  comparing  our  method  with  any  other  fast  wavelet  algorithms  accomplished 
through  convolutions  [6],  which  have  complexity  0(N  log  N),  we  note  that  in  addition  to 
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the  above  mentioned  advantages,  some  processes,  such  as  upsampling  and  downsampling , 
are  eliminated.  Details  of  this  will  be  discussed  in  next  section. 


In  this  section,  we  provide  some  application  examples  of  wavelet  transforms  using 
the  p-product.  The  wavelet  transform  has  become  a very  popular  tool  in  image  processing 
[3,  6,  16,  17].  These  applications  can  be  accomplished  by  use  of  the  fast  algorithm  of 
wavelet  transforms  defined  in  Theorem  6.3.2  with  some  unique  advantages.  Although  we 
provide  only  a few  special  examples,  other  transform  techniques  can  be  accomplished 
in  a similar  fashion. 

6.4.1.  Fourier- Wavelet  Matrix  Expansion 

If  the  characterization  Haar  matrix  of  a wavelet  matrix  a of  rank  m with  genus  g is 
the  finite  Fourier  transform  matrix  fim,  then  the  wavelet  matrix  expansion  for  a function 
/ : Z — > C is  called  a Fourier-wavelet  matrix  expansion. 

If  g= 1,  then  a = fL  In  this  case,  the  Fourier  wavelet  matrix  expansion  for  an 
iV-length  function  f has  the  form 

oo  f m — 1 'j 

f(mr  + l)=  ^ ] 5^  fk(s)wT' * > , r = 0, . . . , m — 1,  (365) 

k—  — oo  v s=0  J 

where  w = e2x,/m.  The  inner  sum  is  the  finite  Fourier  transform  for  the  restriction  of  f 
to  the  interval  ( mk,mk  + m — 1).  That  is  to  say, 


For  g=l,  the  Fourier-wavelet  series  represents  f as  a short  time  Fourier  series;  i.e.,  a 
sequence  of  finite  Fourier  series  for  non-overlapping  windows  of  length  m.  By  the 


6.4.  Applications 


m— 1 


(366) 
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formulae  defined  in  Theorem  6.3.2,  the  Fourier-wavelet  matrix  expansion  under  the 
generalized  matrix  product  is 


where 


f = (f(5)®w*)> 

0 <s<m 


f = 


f(0),f(!),...,f(m  - 1) 


T / 


(367) 


(368) 


?W=  f/o(»)./iW /*-,(*)l  =-(f©mw»), 

L m 772 


and  ws  = 


w°,w\ 


T I 


(369) 


Note  that  both  tensor  product  and  m-product  appearing  in  Equations  (367)  and  (369), 
respectively,  can  be  executed  in  parallel.  This  simplifies  the  computation  of  the  wavelet 
transform  and  its  inverse.  With  sufficient  parallelism,  the  complete  computation  requires 
only  m operations  instead  of  the  2 Nm  operations  used  in  Equations  (365)  and  (366). 


6.4.2.  Wavelet  Decompositions  and  Reconstructions 


In  this  subsection,  we  first  review  algorithms  from  [5]  for  wavelet  decompositions 
and  reconstructions,  and  then  apply  Theorem  6.3.2  to  produce  new  decomposition 
algorithms. 

The  idea  of  a multiresolution  analysis  is  to  write  L2— functions  as  a limit  of 
successive  approximations,  each  of  which  is  a smoothed  version  of  the  original  function. 
The  successive  approximations  thus  construct  different  resolutions,  whence  the  name 
multiresolution  analysis. 

Any  wavelet  generates  a direct  sum  decomposition  of  L2(R),  but  most  of  these 
approximations  are  naturally  associated  with  orthonormal  (o.n.)  wavelet  bases  since 
they  have  good  approximation  properties  [5].  However,  every  o.n.  wavelet  basis  is 
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associated  with  a multiresolution  analysis,  where  L2(R)  is  decomposed  into  a chain  of 
closed  successive  approximation  subspaces 


. . . V_2  C V-i  C Vo  C V\  C V2  . . . 


(370) 


such  that  n jVj  = {0}  and  UJVj  = L2(R).  If  Wj  is  the  orthogonal  complement  of  V3  in 
Vj+ 1,  then  L2(R)  = -j -jWj,  where  denotes  the  direct  sum.  As  pointed  out  in  Section 
1 of  this  chapter,  the  o.n.  wavelet  bases  for  L2(R)  are  formed  from  a function  ip(x), 
such  that  {iftj}k\k  G Z } , where  ipjtk(x)  = V^ipiVx  — k ) spans  Wj.  Associated  with 
the  wavelet  xp  is  another  function  <p,  the  scaling  function,  such  that  {<f>hk\k  € Z},  where 
(f>Jtk(x)  = V^2(f>{V x — k ) forms  a basis  for  Vj.  It  has  been  shown  that  every  function  / 
in  L2(R)  can  be  approximated  as  closely  as  is  desired  by  a function  /jv  G Vjv  for  some 
N G Z [5].  Since  Vj  — Vj-i+Wj-i,  for  any  j G Z,  /jv  has  a unique  decomposition: 


where  f^-i  G VV-i  and  gu~\  G W^-x-  By  repeating  this  process,  we  may  write 


where  fj  G Vj  and  gj  G Wj  for  any  j,  and  M is  chosen  so  that  /n-m  is  sufficiently 
blurred.  The  unique  decomposition  in  the  above  equation  is  called  the  wavelet 
decomposition ; and  the  blur  is  measured  in  terms  of  the  variation  (or  more  precisely, 
frequency  or  number  of  cycles  per  unit  length)  of  Jn-m-  An  easier  stopping  criterion, 
but  less  accurate  than  measuring  variation,  is  to  require  W/n-mW  to  be  smaller  than 
some  threshold. 

In  the  following  paragraphs,  we  will  discuss  an  algorithmic  approach  for  expressing 
fp  as  a direct  sum  of  its  components  gp-i, . . . ,gN-M > and  fN-M  and  recovering 
from  these  components. 


In  = Jn- i + 9N- 1 


(371) 


fN  = 9N- 1 + 9N- 2 + ■ • • + 9N-M  + In~M > 


(372) 
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Let  / 2 denote  the  space  of  all  square-summable  bi-infinite  sequences;  that  is 
{cn}  G / 2 if  and  only  if  ^ \cn\2  < oo.  Since  both  the  scaling  function 

— oo<n<oo 

<j>  G Vo  and  the  wavelet  0 G Wo  are  in  Vi,  and  since  Vi  is  generated  by  <f>i^  = 
21/20(2x  — k ),  where  k G Z,  there  exist  two  sequences  {p*}  and  {^}  £ /2  such  that 

<f>(x)  = y^p*0(2x  - fc),  (373) 

fc 

and 

0(®)  = 9*^(2x  “ *0>  (374) 

j fc 

for  all  a;  G R.  The  above  formulas  are  called  the  two-scale  relations  of  the  scaling 
function  and  wavelet,  respectively.  On  the  other  hand,  since  both  <f>(2x)  and  <f>( 2x  — 1) 
are  in  Vi  and  Vi  = Vo+Wo,  there  are  four  / 2 sequences  which  we  denote  by 
{a_2jt},  {b-2k},  {ai-2fc}>  and  {&i_2*},  with  k G Z,  such  that 

0(2®)  = ^2  [a_2*0(z  — k)  + b_2k^{x  - k)],  and 

k „ (375) 

0( 2x  - 1)  = 2^  [al-2 k<f>(x  - k)  + &1-2 ki’ix  - *)], 
k 

for  all  x G R.  The  above  two  formulas  can  be  combined  into  the  single  formula: 

0(2®  — /)  = ^ [a/-2*0(*  - k)  + &i_2jfc0(z  - k)],  l G Z,  (376) 

k 

which  is  called  the  decomposition  relation  of  0 and  0.  We  now  have  two  pairs  of 
sequences  ({ pk },  {<?*})  and  ({a*},  {&*}),  all  of  which  are  unique  due  to  the  direct  sum 
relationship  Vi  = Vo+Wo.  The  sequences  {pk}  and  {qk}  are  called  reconstruction 
sequences,  while  sequences  {a^}  and  {bk}  are  called  decomposition  sequences.  The 
reason  for  these  names  is  that  the  sequences  can  be  used  to  formulate  reconstruction  and 
decomposition  algorithms. 
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To  describe  these  algorithms,  let  us  first  recall  that  both  fj  e Vj  and  g3  e Wj  have 
unique  series  representations: 


fj(x)  = Y x ~ k) 

k 

with  = jc£  j € l2, 


(377) 


and 


9j{x ) = £<fjW2J'x  ~ *) 
k 

with  dJ  = {4}  £ /2, 


(378) 


where  we  have  intentionally  suppressed  the  normalization  coefficient  2^ 2 by  writing 
out  <f>(2?x  - k)  and  ^»(2Jz  — fc)  instead  of  using  <f>j ^ and  V’y.Jb-  This  allows  us  to  drop 
the  unnecessary  multiple  \/2  in  our  algorithms.  In  the  following  decomposition  and 
reconstruction  algorithms,  the  functions  fj  and  g3  are  represented  by  the  sequences  c; 
and  dJ  as  defined  above. 


Decomposition  Algorithm.  By  applying  Equations  (376)-(378),  we  may  formulate 

Cfc  = Y a/— 2fcc/  1 


(379) 


Observe  that  both  c-7-1  and  d-7-1  are  obtained  from  cJ  by  moving  average  schemes, 
using  the  decomposition  sequences  as  weights , with  the  exception  that  these  moving 
averages  are  sampled  only  at  the  even  integers.  This  is  called  downsampling.  Therefore, 
each  of  the  arrows  in  Figure  6 indicates  a moving  average  followed  by  a downsampling 
at  even  indices.  By  the  auto-stride  property  of  the  /7-product  [24],  a downsampling  is  not 
necessary  in  the  /7-product  wavelet  transform  formula  which  we  will  show  next. 


dN-1  d/V-2  <\N-M 


Figure  6.  Wavelet  decomposition 
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By  taking  m — 2 in  Theorem  6.3.2,  the  wavelet  matrix  a has  the  form 


a = 


do 

bo 


a2g—l 
b2g-l  J 


(380) 


for  some  integer  g.  This  means  that  a*.  = 6*  = 0 unless  0 < k < 2g.  Applying  this 
to  Equation  (379),  we  obtain 


2<7+2fc— 1 

= EaZ-2  ktf  = E a/-2fcc/5 
l l— 2k 


(381) 


2g+2k-l 

= E&Z-2it<4  = E bi_2kCi. 

I l=2k 


Let  c3  — 


Cq,  c[,  . . . , cL  _Js_1  . Then,  by  Theorem  6.3.2,  the  wavelet  decomposition 


2(fl-l) 

algorithm  under  the  2-product  is 

-i.  1 


c;  1 = 


c>0 


ci  ;c^' 

a2g-2 
. a2g—l  _ 


+ • • • + <^-s+2  02 


d'"1  = 

CJ  ©2 

where  cl  = 


dir1^1 


hg-2 

lb2g-l 


+ • ■ ■ + cis+2  ©2 


02 

03 


; cl  ©2 


; cj  ©2 


00 

01 


+ • • ■ + <4-1  ©2 


a2g-2 
.02,7-1  J. 


(382) 


+ • • • + ca-l  ©2 


hg-2 

b2g-\ 


is  of  the  same  size  as  c3 . 


-2  V °2i+l> 

Note  that  if  the  length  of  c3  is  N,  then  computing  cJ_1  by  using  Equation  (379) 
requires  2g-  (y  + g — l)  = Ng  + 2g(g  - 1)  operations.  But  if  we  employ  Equation  (382) 
instead,  then  each  2-product  needs  2(g-l)  operations  for  computing  elements  of 
while  each  2-product  needs  N operations  for  computing  elements  of  c+_1.  Therefore,  the 
maximum  execution  time  for  complete  parallel  processing  is  no  more  than  N operations, 
which  means  that  it  requires  at  least  g times  fewer  computations.  Additionally,  in  real 
implementation,  by  employing  the  advantages  of  the  variations  of  the  length  of  the  support 
of  each  cl,  more  time  can  be  saved.  As  pointed  out  in  Section  3,  if  there  are  enough 
processors  on  parallel  machines,  the  execution  time  is  2 g vs  Ng  -f  2g(g  — 1). 
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Reconstruction  Algorithm.  By  applying  Equations  (373),  (374),  (377)  and  (378),  we 
have 


4 = \pk-21cj  1 + qk-2i4  1 


(383) 


Here,  cJ  is  obtained  from  c-7-1  and  dJ_1  by  two  moving  averages,  using  the 
reconstruction  sequences  as  weights,  with  the  exception  that  an  upsampling  is  required 
before  the  moving  averages  are  performed.  More  precisely,  in  Equation  (383),  the  indices 
of  the  input  sequences  {c^— 1 } and  {dj-1}  316  multiplied  by  2,  and  zeros  are  used  as  the 
terms  with  odd  indices  in  the  new  input  sequences,  where  the  (discrete)  convolutions 
are  taken  with  respect  to  {pm}  and  {qm}-  Therefore,  each  of  the  arrows  in  Figure  7 
indicates  a moving  average  followed  by  upsampling.  However,  in  what  follows,  we 
will  show  that  the  upsampling  is  not  necessary  when  one  applies  the  p-product  wavelet 
transform  formula  and  its  unique  properties  [24]. 


N-M 

d 

N-M 


A/-/W+1 

d 

N-M*!  \ 

c ► 


N 

c 


f 


Figure  7.  Wavelet  reconstruction 


Similar  to  the  decomposition  algorithm,  by  taking  m- 2,  the  wavelet  matrix  a has 
the  form 


a = 


P o 

qo 


P2g-l 
q2g-\  \ ’ 


(384) 
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for  some  integer  g.  Thus,  by  Theorem  6.3.2,  the  wavelet  reconstruction  algorithm  under 
the  p-product  is 

cj  = C*-1  ® [po  Pi]  + Cj_1  ® [P2  P3]  + ■ • . + cj:}  <g>  \p2g-2  P2g-l] 

+dJ_1  ® [90  9l]  + dj-1  ® [q2  ©]  + •••  + djl}  ® [ q2g-2  92S-l], 
where  c;--1  = , ci~+1, . . . is  of  the  same  size  as  cJ-1, 


and  d;  1 = 


dJ~l  d]~l 

a-i  >“-,+!>••• 


is  of  the  same  size  as  d-7 


-1 


(385) 


6.4.3.  Pyramid  Transforms 

The  pyramid  transform  is  a special  case  of  multiresolution  analysis,  which 
employs  non-overlapping  window  functions  to  transform  an  image  from  one  resolu- 
tion to  another.  An  M-pyramid  (or  pyramid  in  the  matrix- sequence  sense)  P is  a 
sequence  {M(L),  M(L  - 1), . . . , Af(O)}  of  arrays  where  M(L ) represents  an  input 
image,  M(i  — 1)  is  a version  of  M(i)  at  half  the  resolution  of  M(i)  and  so  on, 
with  M(0)  a single  pixel.  That  is  if  M(L ) is  2L  x 2L,  the  successive  layers  are 

2l~ 1 x 2l~1,2l~2 * * * *  x 2l~2,...,2  x 2,1  x 1.  So  that  the  height  of  the  pyramid  is  L. 

A pixel  of  an  image  M(i ) is  referred  to  as  a cell  of  the  pyramid  at  level  i.  Each  cell 

has  (in  general)  four  brothers  in  its  own  layer,  four  children  in  the  layer  below  it,  and 

one  parent  in  the  layer  above.  Arrays  to  be  accepted  are  input  to  the  base,  and  the  apex 

cell  is  the  distinguished  cell. 

The  pyramid  structure  can  be  used  for  image  search,  edge  detection,  pattern 
recognition,  etc.  The  simplest  application  of  pyramids  is  to  provide  reduced-resolution 
versions  of  an  image.  The  processing  time  to  transform  a 64  x 64  array  is  one-fourth  the 
time  for  a 128  x 128  array  in  most  circumstances.  When  several  different  operations  are 
to  be  performed  on  an  image,  each  operation  should  use  only  the  amount  of  resolution 
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necessary.  A pyramid  provides  a representative  selection  of  the  possible  versions  of  an 
image  at  different  degrees  of  resolution. 

Since  the  pyramid  transform  uses  only  non-overlapping  window  functions,  it  is 
advantageous  to  employ  the  p-product.  In  what  follows,  we  show  how  the  p-product 
allows  us  to  traverse  a pyramid  in  both  directions. 

Suppose  we  have  a pyramid  with  base  size  2"  x 2";  that  is  the  input  image 

I'M 


a £ F2"x2n  and  b — 


£ F2xi-  Then 


b2 


C2  — [(a©2^)  ©2^]  £ F2n-lx2n-l, 


(386) 


with 


C2 ij  = (a(2i-l)(2j-l)  obl°  ^l)7(«(2.-l)(2j)7bl  o b2)  7 
(a(2.)(2i-l)  0 0 ^2)7(a(2»)(2i)  0 0 h)- 


Therefore, 


and 


c4  = 


(((a@2b)'©2b)'®2b)  ®2b 

(...(a®2b)'®2...)'®2b]' 


£ F2»-2x2”-2J 


2n~ix2n~i ' 


It  follows  that 


C 2n=  (•  ..(a©2b),©2---/®2b] 


£ F. 


(387) 


(388) 


(389) 


(390) 


In  particular,  if  7=V,  then  C2n  = Va  = max  (a),  where  max  (a)  denotes  the 
maximum  value  of  a.  Similarly,  if  7 = A,  then  C2n  = Aa  = min(a)  by  choosing 


b = 

C2n  : 


. Choosing  the  operation  @p  instead  of  El  P and  b = 

2^Ea- 


1/2 

L1/2J 


, we  obtain 
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In  general,  if  we  have  a pyramid  with  base  size  kn  x kn  and 

h 


b = 


bk 


€ Fjfcxi, 


then  for  a £ Fk”xkn>  we  obtain 


ckn=  (...(a^b)'©*...)'©^] 


£ F. 


(391) 


(392) 


The  above  discussion  shows  how  to  traverse  a pyramid  from  bottom  to  top.  We 
now  present  scheme  for  traversal  from  top  to  bottom. 

£ F 2x2?  and  p=l.  Then 


b\\  b\2 
&21  &22 


Let  a = [ai]  £ F,  b = 

Ci  = a ® b = 


a\  o fen  a i o fei2 
a\  o &2i  a\  o &22 


£ F 


2x2- 


and 


c2  = 

(a  ® b)  (g>  b 

" (ai  o feia)  o fen 

(ai  o fen)  o fe12 

(oi  o fe12)  o fen 

(d  ofe12)  ofe12* 

(oi  o fen)  o fe21 

(ai  o fexl ) o fe22 

(ai  o fe12)  o fe21 

(ai  o fe12)  o fe22 

(ai  o fe2i)  o fen 

(ai  o fe2i)  o fe12 

(ai  o fe22)  o fen 

(ai  o fe22)  o fei2 

_(ai  o fe21)  o fe21 

(ai  o fe21)  o fe22 

(ai  o fe22)  o fe21 

(ai  o fe22)  o fe22 . 

Thus 


(393) 


£ F 


4x4" 


(394) 


c„  = (. . . (a  <g>  b)  ® b . . . ® b)  £ F2«x2»,  for  some  n £ Z+. 


(395) 


If  we  take  b = 


If 


1 1 
1 1 


, and  o=multiplication,  then  if 
a = [ai],  ci  = 


a i ai 
ai  a,\ 


a = 


an  an 
021  «22 


£ F2x2 


(396) 


(397) 
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then 


Cl  = 


Hence,  if  a G F2;x2;,  then 


a it 

Oil 

012 

012 

an 

Oil 

«12 

012 

021 

021 

022 

022 

«21 

021 

022 

022 

(398) 


cn  = (•  • • (a  ® b)  <g)  b . . . ® b)  G F2;+r»x2<+n. 


(399) 


In  general  if  a G F/Xm,  bGF„x?,  then  for  any  positive  integer  i we  obtain 


c.‘  = (•  • • (a  ® b)  ® b . . . ® b)  € F(/xn,)x(rn><gi). 


(400) 


CHAPTER  7 

CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


We  have  defined  a new  matrix  product — the  generalized  matrix  product  (GMP) 
or  /7-product.  We  demonstrated  that  this  new  operation  possesses  a rich  mathematical 
content  and  is  a useful  tool  for  the  development  of  parallel  algorithms  for  a large  class 
of  orthogonal  and  nonorthogonal  transforms.  It  can  be  used  to  express  a wide  variety 
of  operations  in  a common  framework. 

More  specifically,  we  have  shown  that: 

(1)  Many  properties  of  matrix  product  and  tensor  product  theory,  such  as  the 
associative  law,  distributive  law,  determinant  theory,  etc.,  have  their  counterparts 
in  the  /7-product  theory.  In  addition,  there  are  some  characteristics  inherent  in  the 
p-product,  such  as  the  parallel  property  (Theorem  2.2.10),  Theorem  2.2.8,  etc., 
that  make  the  /7-product  a useful  tool  in  parallel  computing. 

(2)  The  linear  algebra  generalized  by  the  GMP  provides  a rich  mathematical  back- 
ground for  this  new  operation.  Just  as  for  linear  transform  representation, 
the  matrix  representation  of  the  /7-product  transformation  permits  us  to  employ 
matrix  properties  to  study  /7-product  transformations.  The  applications  of  these 
transformations  in  signal  processing  remain  to  be  explored. 

(3)  The  close  relationship  between  the  tensor  product  and  the  /7-product  allows  us 
to  develop  applications  of  the  /7-product  in  particle  physics  and  matrix  calculus, 
as  well  as  in  various  fields  of  matrix  theory  in  which  the  tensor  product  plays  a 
significant  role  [13]. 
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(4)  The  p-product  furnishes  an  alternative  and  extensive  algebraic  formulation  of 
a class  of  linear  transformations.  It  can  be  applied  to  express  various  image 
processing  transforms.  It  provides  efficient  computations  for  such  well  known 
transforms  as  the  Fourier,  Walsh  and  generalized  Walsh  transforms,  as  well  as 
simplifications  of  algorithms  for  their  computation  due  to  replacement  of  large 
blocks  of  code  by  short  algebraic  p-product  expressions.  We  have  also  shown  that 
it  is  possible  to  use  the  p-product  to  establish  simplified  algorithms  for  efficient 
matrix  operations  over  a class  of  generalized  tensor  matrices.  Thus,  any  orthogonal 
or  nonorthogonal  transforms  obtainable  from  a series  of  matrix  tensor  products 
can  be  expressed  concisely  in  terms  of  the  p-product,  and  may  be  efficiently 
implemented  for  computer  applications. 

(5)  The  generalized  matrix  product  can  successfully  decompose  a long  summation  into 
short  ones,  yielding  efficient  computation  in  parallel  machines.  It  lends  itself  well 
to  expressing  the  wavelet  transform  in  matrix  form  and  provides  simple  and  fast 
algorithms  for  the  wavelet  transform  and  its  inverse  using  the  parallelism  inherent 
in  the  p-product  formulation. 

The  generalized  matrix  product  is  a very  recent  discovery.  Thus,  there  is  plenty 
of  room  for  additional  research.  We  now  list  some  suggestions  for  further  research.  We 
have  attempted  to  present  the  problems  in  the  same  order  as  the  material  in  the  main 
body  of  the  dissertation  on  which  they  are  based. 

Even  though  certain  properties  of  the  p-product  have  been  developed  in  this 
dissertation,  they  were  based  only  on  its  utilities  related  to  the  research  presented  here. 
A wealth  of  mathematical  results,  with  possible  applications  to  real  world  problems,  are 
obviously  obtainable  from  investigating  other  properties  of  this  product.  The  successful 
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application  of  the  parallel  property  (Theorem  2.2.10)  in  wavelet  transform  presented  in 
this  document  is  a clear  example. 

We  have  mapped  the  matrix  representation  theorem  to  the  p-product.  However,  in 
the  theory  of  linear  transformations,  each  homomorphism  can  be  represented  by  different 
matrices  relative  to  different  bases,  does  the  same  property  hold  for  the  theory  of  p-product 
transformations?  How  about  the  other  linear  algebra  theorems?  Can  we  also  map  them  to 
this  product?  Are  there  some  p-product  transforms  which  are  as  important  as  regular  linear 
transforms?  What  are  their  uses?  The  generalized  vector  space  product  we  presented 
here  is  similar  to  the  tensor  product  vector  space.  What  are  its  applications?  Can  other 
techniques  in  linear  algebra,  such  as  eigenvalue/eigenvector,  SVD,  rank,  solutions  to 
systems  of  equations,  etc.  be  described  in  a similar  way  using  the  p-product? 

In  Chapter  4 and  5,  we  have  applied  the  p-product  into  Fourier  and  Walsh 
transforms.  The  p-product  presentations  of  the  generalized  transforms , which  contain 
both  Fourier  and  Walsh  transforms  as  special  cases,  have  yet  to  be  established. 

The  last  question  concerns  the  applications  of  the  p-product  in  wavelet  transform 
theory.  We  have  established  a fast  wavelet  transform  in  terms  of  the  p-product  based  on 
a given  wavelet  matrix.  Is  this  possible  to  apply  the  p-product  in  deriving  mathematical 
theorems  of  wavelet  systems?  How  successful  will  it  be  if  we  apply  the  p-product  wavelet 
transform  to  solve  the  differential  equations?  All  these  applications  may  involve  the 
development  of  the  p- product  properties.  In  addition,  recall  that  the  efficient  computation 
of  the  fast  wavelet  transform  algorithm  presented  here  relies  heavily  on  implementations 
and  parallel  machines,  finding  out  an  optimal  result  is  important. 

In  the  original  paper  of  the  generalized  matrix  product  [23],  Ritter  pointed  out  that 
this  new  product  forms  a basis  for  the  structure  of  image  algebra.  What  is  the  connection 
between  these  two  topics? 


APPENDIX 

COMPUTER  PROGRAM 


In  this  appendix,  we  present  a computer  program  used  to  implement  p -product  fast 
Fourier  transform  under  Matlab  3.0  for  computing  Fourier  transform  of  a vector  f with  the 
length  being  the  powers  of  2.  It  is  obtained  by  modifying  the  algorithm  (CTFFT)  based 
on  Cooley-Tukey  method  and  provided  by  Prof.  Sigmon.  The  function  BRACEWELL 
is  a subroutine  of  the  function  CTFFT,  which  performs  the  reordering  operation.  Two 
functions  (PFFT  and  CTFFT)  are  almost  same  except  PFFT  employs  the  GMP  idea, 
which  requires  no  bit-reversal  permutation  (function  BRACEWELL).  As  the  result  of 
this,  PFFT  not  only  makes  the  program  shorter  but  also  saves  up  to  30%  of  the  CPU 
time  as  shown  on  Table  1. 
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function  a = PFFT(f) 

a = f(:); 

n = length(f); 
n2  = n/2; 
w — exp  (— 7rz'/n2); 
u = ones(n2, 1); 
for  k = 2 : n2, 

u(k ) = u(A:  — 1)  * w, 

end 


bw  = 1;  k = n; 
while  k > 1 

k2  — k/ 2;  bw2  = bw  * 2; 
bw3  = bw  + 1;  Ar3  = A:2  -+- 1; 
c = a(l  : k2 , 1 : bw ); 
c(l  : k2,bw3  : bw2)  = c; 
d = a(&3  : k,  1 : bw)- 
d(l  : k2,bw3  : bw 2)  = — d; 
a = c -|-  d; 

if  k > 3 

w = u(l  : bw  : n2); 
w = w(:,  ones(l,  bw))-, 
a(:,  bw3  : bw2)  = a(:,  bw3  : bw2).  * w; 
end 

bw  = bw2 ; k = k2\ 

c = 0;  d = 0; 


end 
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function  a = CTFFT(f) 


a = f(0; 

n = length(f);  n2  = n/ 2; 
w = exp  (— ni/n2); 
u = ones(n2, 1); 
for  k = 2 : n2, 

u(k)  = u (k  — 1)  * w, 

end 

a = bracewell(a.)] 

bw  = 1;  k = n; 
while  k > 1 
k2  = k/2; 

w = u(l  : k2  : n2); 
w = w(:,ones(l,  k2)); 
c = 2ero.s(bw,  k); 
c(:)  = a;  b = c; 
c(:,2  : 2 : k)  = c(:,l  : 2 : k); 
b(:,  2 : 2 : k)  = — w.  * b(:,  2 : 2 : k); 
b(:,l  : 2 : k)  = -b(:,2  : 2 : k); 
a = c + b;  a = a(:); 
bw  = 2 * b2;  k = k2; 


end 
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function  a = bracewell(a) 

n = length(  a); 
m = round(log(n)/log(  2)); 
sqrtn  = <Ifloor{{m  + 1) /2); 

indices  = 0 : sqrtn  — 1; 
b = 2; 

while  b < sqrtn 

twice  = 2 * indices^  1 : 6); 
b — 2 * b; 

indices^  1 : 6)  = [twice,  twice  + 1]; 

end 

if  rem(m,  2)  ==  1,  b = 6/2;  end 

for  q = 2 : 6 

r — b*  indices(q)  + 1; 
np  = q : b : r; 

rp  = r + indices^  1 : indices(q))\ 
a ([np  rp])  = a([rp  np]); 


end 
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